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PREFACE 


This  volume  is  part  of  a  four-volume  set  that  summarizes  the  research  of  participants  in 
the  1999  AFOSR  Summer  Research  Extension  Program  (SREP).  The  current  volume. 
Volume  1  of  4,  presents  the  final  reports  of  SREP  participants  at  Armstrong  Laboratory. 

Reports  presented  in  this  volume  are  arranged  alphabetically  by  author  and  are  numbered 
consecutively  -  e.g.,  1-1,  1-2,  1-3;  2-1,  2-2,  2-3,  with  each  series  of  reports  preceded  by 
a  35  page  management  summary.  Reports  in  the  four-volume  set  are  organized  as 
follows: 


VOLUME  TITLE 

1  Armstrong  Research  Laboratory 

2  Phillips  Research  Laboratory 

3  Rome  Research  Laboratory 

4  Wright  Research  Laboratory 
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1999  SUMMER  RESEARCH  EXTENSION  PROGRAM  (SREP)  MANAGEMENT  REPORT 


1.0  BACKGROUND 


Under  the  provisions  of  Air  Force  Office  of  Scientific  Research  (AFOSR)  contract  F49620-90-C- 
0076,  September  1990,  Research  &  Development  Laboratories  (RDL),  an  8(a)  contractor  in 
Culver  City,  CA,  manages  AFOSR’s  Summer  Research  Program.  This  report  is  issued  in  partial 
fulfillment  of  that  contract  (CLIN  0003 AC). 

The  Summer  Research  Extension  Program  (SREP)  is  one  of  four  programs  AFOSR  manages 
under  the  Summer  Research  Program.  The  Summer  Faculty  Research  Program  (SFRP)  and  the 
Graduate  Student  Research  Program  (GSRP)  place  college-level  research  associates  in  Air  Force 
research  laboratories  around  the  United  States  for  8  to  12  weeks  of  research  with  Air  Force 
scientists.  The  High  School  Apprenticeship  Program  (HSAP)  is  the  fourth  element  of  the  Summer 
Research  Program,  allowing  promising  mathematics  and  science  students  to  spend  two  months  of 
their  summer  vacations  working  at  Air  Force  laboratories  within  commuting  distance  from  their 
homes. 

SFRP  associates  and  exceptional  GSRP  associates  are  encouraged,  at  the  end  of  their  summer 
tours,  to  write  proposals  to  extend  their  summer  research  during  the  following  calendar  year  at 
their  home  institutions.  AFOSR  provides  funds  adequate  to  pay  for  SREP  subcontracts.  In 
addition,  AFOSR  has  traditionally  provided  further  funding,  when  available,  to  pay  for  additional 
SREP  proposals,  including  those  submitted  by  associates  from  Historically  Black  Colleges  and 
Universities  (HBCUs)  and  Minority  Institutions  (Mis).  Finally,  laboratories  may  transfer  internal 
funds  to  AFOSR  to  fund  additional  SREPs.  Ultimately  the  laboratories  inform  RDL  of  their 
SREP  choices,  RDL  gets  AFOSR  approval,  and  RDL  forwards  a  subcontract  to  the  institution 
where  the  SREP  associate  is  employed.  The  subcontract  (see  Appendix  1  for  a  sample)  cites  the 
SREP  associate  as  the  principal  investigator  and  requires  submission  of  a  report  at  the  end  of  the 
subcontract  period. 

Institutions  are  encouraged  to  share  costs  of  the  SREP  research,  and  many  do  so.  The  most 
common  cost-sharing  arrangement  is  reduction  in  the  overhead,  fringes,  or  administrative  charges 
institutions  would  normally  add  on  to  the  principal  investigator’s  or  research  associate’s  labor. 
Some  institutions  also  provide  other  support  (e.g.,  computer  run  time,  administrative  assistance, 
facilities  and  equipment  or  research  assistants)  at  reduced  or  no  cost. 

When  RDL  receives  the  signed  subcontract,  we  fund  the  effort  initially  by  providing  90%  of  the 
subcontract  amount  to  the  institution  (normally  $18,000  for  a  $20,000  SREP).  When  we  receive 
the  end-of-research  report,  we  evaluate  it  administratively  and  send  a  copy  to  the  laboratory  for  a 
technical  evaluation.  When  the  laboratory  notifies  us  the  SREP  report  is  acceptable,  we  release 
the  remaining  funds  to  the  institution. 


Introduction  -1 


2.0  THE  1999  SREP  PROGRAM 


SELECTION  DATA:  A  total  of  381  faculty  members  (SFRP  Associates)  and  130  graduate 
students  (GSRP  associates)  applied  to  participate  in  the  1998  Summer  Research  Program.  From 
these  applicants  85  SFRPs  and  40  GSRPs  were  selected.  The  education  level  of  those  selected 
was  as  follows: 


1998  SRP  Associates,  by  Degree 

SFRP 

GSRP 

PHD 

MS 

MS 

BS 

83 

1 

3 

19 

Of  the  participants  in  the  1998  Summer  Research  Program  65  percent  of  SFRPs  and  20  percent 
of  GSRPs  submitted  proposals  for  the  SREP.  Fifty-four  proposals  from  SFRPs  and  eleven  from 
GSRPs  were  selected  for  funding,  which  equates  to  a  selection  rate  of  65  %  of  the  SFRP  proposals 
and  of  20%  for  GSRP  proposals. 


1999  SREP:  Proposals  Submitted  vs.  Proposals  Selected 

Summer 

1998 

Participants 

Submitted 

SREP 

Proposals 

SREPs 

Funded 

SFRP 

85 

54 

34 

GSRP 

40 

11 

2 

TOTAL 

125 

65 

36 

The  funding  was  provided  as  follows: 


Contractual  slots  funded  by  AFOSR  36 

Laboratory  funded  _0 

Total  36 

Four  HBCU/MI  associates  from  the  1998  summer  program  submitted  SREP  proposals;  four  were 
selected  (none  were  lab-funded;  all  were  funded  by  additional  AFOSR  funds). 
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Proposals  Submitted  and  Selected,  by  Laboratory 

Applied 

Selected 

Armstrong  Research  Site 

15 

3 

Air  Logistic  Centers 

1 

0 

Arnold  Engineering  Development  Center 

1 

0 

Phillips  Research  Site 

10 

8 

Rome  Research  Site 

12 

7 

Wilford  Hall  Medical  Center 

1 

0  1 

Wright  Research  Site 

25 

18 

TOTAL 

65 

36 

The  212  1998  Summer  Research  Program  participants  represented  60  institutions. 


Institutions  Represented  on  the  1998  SRP  and  1999  SREP 

Number  of  schools 
Represented  in  the 
Summer  98  Program 

Number  of  schools 
represented  in 
submitted  proposals 

Number  of  schools 
represented  in 
Funded  Proposals 

60 

36 

29 

Thirty  schools  had  more  than  one  participant  submitting  proposals. 
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The  selection  rate  for  the  60  schools  submitting  1  proposal  (68%)  was  better  than  those 
submitting  2  proposals  (61%),  3  proposals  (50%),  4  proposals  (0%)  or  5+  proposals  (25%). 
The  4  schools  that  submitted  5+  proposals  accounted  for  30  (15%)  of  the  65  proposals 
submitted. 

Of  the  65  proposals  submitted,  35  offered  institution  cost  sharing.  Of  the  funded  proposals 
which  offered  cost  sharing,  the  minimum  cost  share  was  $1274.00,  the  maximum  was 
$38,000.00  with  an  average  cost  share  of  $12,307.86. 


Proposals  and  Institution  Cost  Sharing 

Proposals 

Submitted 

Proposals 

Funded 

With  cost  sharing 

35 

30 

Without  cost  sharing 

30 

6 

Total 

65 

36 

The  SREP  participants  were  residents  of  31  different  states.  Number  of  states  represented  at 
each  laboratory  were: 


States  Represented,  by  Proposals  Submitted/Selected  per  Laboratory 

Proposals 

Submitted 

Proposals 

Funded 

Armstrong  Research  Laboratory 

15 

3 

Air  Logistic  Centers 

1 

0 

Arnold  Engineering  Development  Center 

1 

0 

Phillips  Research  Laboratory 

10 

8  i 

Rome  Research  Laboratory 

12 

7 

Wilford  Hall  Medical  Center 

1 

0 

Wright  Research  Laboratory 

25 

18 

Six  of  the  1999  SREP  Principal  Investigators  also  participated  in  the  1998  SREP. 


ADMINISTRATIVE  EVALUATION:  The  administrative  quality  of  the  SREP  associates’  final 
reports  was  satisfactory.  Most  complied  with  the  formatting  and  other  instructions  provided  to 
them  by  RDL.  Thirty-six  final  reports  have  been  received  and  are  included  in  this  report.  The 
subcontracts  were  funded  by  $897,309.00  of  Air  Force  money.  Institution  cost  sharing  totaled 
$356,928.00. 
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TECHNICAL  EVALUATION:  The  form  used  for  the  technical  evaluation  is  provided  as 
Appendix  2.  Thirty-two  evaluation  reports  were  received.  Participants  by  laboratory  versus 
evaluations  submitted  is  shown  below: 


Participants 

Evaluations 

Percent 

Armstrong  Laboratory 

3 

2 

95.2 

Phillips  Laboratory 

8 

8 

100 

Rome  Laboratory 

7 

7 

100 

Wright  Laboratory 

18 

15 

91.9 

Total 

36 

32 

95.0 

Notes: 

1:  Research  on  four  of  the  final  reports  was  incomplete  as  of  press  time  so  there  aren’t  any  technical 
evaluations  on  them  to  process,  yet.  Percent  complete  is  based  upon  20/21  =95.2% 

PROGRAM  EVALUATION:  Each  laboratory  focal  point  evaluated  ten  areas  (see  Appendix 
2)  with  a  rating  from  one  (lowest)  to  five  (highest).  The  distribution  of  ratings  was  as  follows: 


Rating 

Not  Rated 

1 

2 

3 

4 

5 

#  Responses 

7 

1 

7 

62  (6%) 

226  (25%) 

617  (67%) 

The  8  low  ratings  (one  1  and  seven  2’s  )  were  for  question  5  (one  2)  “The  USAF  should 
continue  to  pursue  the  research  in  this  SREP  report”  and  question  10  (one  1  and  six  2’s)  “The 
one-year  period  for  complete  SREP  research  is  about  right”,  in  addition  over  30%  of  the 
threes  (20  of  62)  were  for  question  ten.  The  average  rating  by  question  was: 


Question 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

Average 

4.6 

mai 

4.7 

wsm 

4.6 

■SB 

4.8 

4.5 

4.6 

4.0 
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The  distribution  of  the  averages  was: 


AREA  AVERAGES 


Area  10  “the  one-year  period  for  complete  SREP  research  is  about  right”  had  the  lowest 
average  rating  (4.1).  The  overall  average  across  all  factors  was  4.6  with  a  small  sample 
standard  deviation  of  0.2.  The  average  rating  for  area  10  (4. 1)  is  approximately  three  sigma 
lower  than  the  overall  average  (4.6)  indicating  that  a  significant  number  of  the  evaluators  feel 
that  a  period  of  other  than  one  year  should  be  available  for  complete  SREP  research. 
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The  average  ratings  ranged  from  3.4  to  5.0.  The  overall  average  for  those  reports  that  were 
evaluated  was  4.6.  Since  the  distribution  of  the  ratings  is  not  a  normal  distribution  the  average 
of  4.6  is  misleading.  In  fact  over  half  of  the  reports  received  an  average  rating  of  4.8  or 
higher.  The  distribution  of  the  average  report  ratings  is  as  shown: 


AVERAGE  RATINGS 


It  is  clear  from  the  high  ratings  that  the  laboratories  place  a  high  value  on  AFOSR’s  Summer 
Research  Extension  Programs. 
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3.0  SUBCONTRACTS  SUMMARY 


Table  1  provides  a  summary  of  the  SREP  subcontracts.  The  individual  reports  are  published 
volumes  as  shown: 


Laboratory  Volume 

Armstrong  Research  Laboratory  1 

Phillips  Research  Laboratory  2 

Rome  Research  Laboratory  3 

Wright  Research  Laboratory  4 
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SREP  SUB-CONTRACT  DATA 


Report  Author  Author’s  Sponsoring  .  Contract  Univ.  Cost 

Author’s  University _  Degree  Lab  Performance  Period  Amount  Share 


Graetz ,  Kenneth  PhD 

Department  of  Psycho logy  99-0803 

University  of  Dayton,  Dayton,  OH 

Kannan ,  Nandini  PhD 

Statistics  99-08  04 


Univ  of  Texas  at  San  Antonio,  San  Antonio,  TX 


Ramesh  ,  Ramaswamy  PhD 

Magemeut  Science/Systems  99-0802 

Research  Foundation  of  SUNY,  Buffalo,  NY 

Le ,  Vanessa  BS 

Biochemistry  99-0805 

Univ  of  Texas  at  Austin,  Austin,  TX 

Gill ,  Gunnm  PhD 

EE  99-0834 

Naval  Postgraduate  School,  Monterey,  CA 

Hinde ,  Robert  PhD 

Physical  Chemistry  99-0801 

Univ  of  Tennessee,  Knoxville,  TN 

Jeffs ,  Brian  PhD 

Electrical  Engineering  99-0828 


Brigham  Young  University,  Provo,  UT 


Leo  ,  Donald  PhD 

Mechanical  &  Aerospace  99-0035 

Virginia  Tech,  Blacksburg,  VA 

Lodhi ,  M.  Arfin  PhD 

Nuclear  Physics  99-0832 

Texas  Tech  University,  Lubbock,  TX 

McHugh.  John  PhD 

Applied  Mechanics  99-0833 

University  of  New  Hampshire,  Durham,  NH 

Steinberg ,  Stanly  PhD 

Mathematics  99-082  9 

University  of  New  Mexico,  Albuquerque,  NM 

Stephens  D  ,  Kenneth  MA 

99-0830 

University  of  North  Texas,  Denton,  TX 

An  as  ,  E rcument  PhD 

Electrical  Engineering  99-0808 

Syracuse  University,  Syracuse,  NY 

Go  pal  an  .  Kaliappan  PhD 

Electrical  Engineering  9  9-0814 


Purdue  Research  Foundation,  West  Lafayette,  IN 


AIVCF  01/01/99  12/31/99  S24983.00  $0.00 

Conflict  resolution  in  distributed  Meetmgl; 
Using  Collaboration  Technology  to 

ALCF  01/01/98  12/31/98  $22492.00  $4478.00 

Altitude  decompression  Sickness:  Modeling  and 
Prediction 

AL/CF  01/01/99  12/31/99  $24979.00  $0.00 

Modeling  and  Analsis  of  DMT  Systems:  Training 
Effectiveness,  Costs,  Resource  Ma 

AL/CF  01/01/98  12/31/98  $25000.00  $0.00 

a  Study  on  Stress  -  Induced  Alterations  In 
Blood-Brain  Barrier  Permeability  to  ?yr 

PL/VT  01/01/99  12/31/99  $25000.00  $0.00 

Adaptive  signal  Processing  and  its  Applications 
in  Space  based  Radar 

WL/PO  01/01/98  12/31/98  $25000.00  $3976.00 

Dopant  - Induced  Infrafed  Activity  in  Sclid 
Hydrogen  An  AB  Initio  and  Quantum  Mont 

PL/LI  01/01/98  12/31/98  $25000.00  $19177.00 

Algebraic  Methods  for  Improved  blimd  Restoration 
of  Adaptive  Optics  Images  of  Sp 

PL/VT  01/01/99  12/31/99  $24999.00  $7416.00 

Self-Sensing  Acoustic  Sources  For  Interior  Noise 
Control  in  Payload  Fairings 

PL/VT  01/01/99  12/31/99  $25000.00  $0.00 

Invest igat ioin  into  Time-Dependent  Power  Losses 
from  AMTEC  Components 

PLVT  01/01/99  12/13/99  $25000.00  $7000.00 

Atmospheric  Gravity  Waves  Near  The  Tropccause 

PL/LI  01/01/99  12/31/99  $25000.00  $0.00 

Lie-Algebraic  Representations  of  Product 
Integrals  of  Variable  Matrices 

PLU  01/01/99  12/31/99  $25000.00  $16764.00 

Simulation  of  a  Magnetized  Target  Fusion  loncept 
Using  MACH2 

WL/AA  01/01/99  12/31/99  $25000.00  $13000.00 

Realization  of  Low  Noise  MMIC  Amplifier  as  a 
Microwave-to-Optics  Link  for  Radar 

RLTR  01/01/99  12/31/99  $25000.00  S38168.00 

Detection  of  Acoustic  Correlates  of  Stress  from 
Modulation  Characteristics 
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SREP  SUB-CONTRACT  DATA 


Report  Author  Author’s 

Author’s  Univenity  Degree 

Hung ,  Donald  PhD 

Electrical  Engineering  99-0812  An 

Washington  State  University,  Richland,  WA 


Sponsoring 

Lab 

RL/IR 


Performance  Period 
1)1/01/99  12/31/99 


Contract 

Amount 

125000.00 


Univ.  Cost 
Share 

$23008.00 


Investigation  on  Accelerating  the  Ray-Tracing 
Comppurations 


Lutoborski ,  Acton 

PhD 

Applied  Mathematics 

99-0811 

Syracuse  Univenity,  Syracuse,  NY 

Panda ,  Brajendra 

PhD 

Computer  Science 

99-0810 

University  of  North  Dakota,  Grand 

Forks,  ND 

RLTO  01/01/99  12/31/99  $25000.00  $1274.00 

Transform  Methods  for  Watermarking  Digital 
Images 

RL/TR  01/01/98  12/31/98  $24942.00  $2600.00 

Implementation  of  Petri  Sets  Based  Multi-source 
Attack:  Detection  Model 


Potter ,  Jerry  PhD 

Computer  Science  99-08  09 

Kent  State  Univenity,  Kent,  OH 

Upadhyaya ,  Shuabhu  PhD 

Elec  &  Comp  Engineering  99-0813 

SUNY  Buffalo,  Buffalo,  NY 

Ahmed ,  Farid  PhD 

Electrical  engineering  99-0806 

Penn  State  Uni-Erie,  Erie,  PA 

Belfield ,  Kevin  PhD 

Chemistry  99-0816 

University  of  Central  Florida,  Orlando,  FL 

Buck ,  Gregory  PhD 

Mechanical  Engineering  99-0818 

S  Dakota  School  of  Mines/Tech,  Rapid  City,  SD 

Gilcrease ,  Patrick  PhD 

Chemical  Engineering  99-0815 

University  of  Wyoming,  Laramie,  WY 

Johnson ,  Jeffrey  PhD 

Electrical  Engineering  and  99-0823 

University  of  Toledo,  Toledo,  OH 

Kapila ,  Vikram  PhD 

Aerospace  engineering  99-0820 

Polytechnic  Inst  of  New  York,  Brooklyn,  NY 

Kihm  ,  Kenneth  PhD 

Mechanical  Engineering  99-0821 

Texas  Engineering  Experiment  Station,  College 

Li  ,  Rongxing  PhD 


Photogrammertry  &  Remote  Sensing  99-0831 
Ohio  State  Univenity  ,  Columbus,  OH 


Lin  ,  Chun -Shin  PhD 

Electrical  Engineering  99-0826 

Univ  of  Missouri  -  Columbia,  Columbia,  MO 

Liu  ,  Chaoqun  PhD 

Applied  Mathematics  99-0819 


Louisiana  Tech  University,  Ruston,  LA 


RL/IR  01/01/99  12/31/99  $25000.00  $52767.00 

Algorithms  for  Data  Intensive  Knowledge 
Discovery 

RL/IR  01/01/98  12/31/98  $25000.00  $6430.00 

a  Distrubuted  concurrent  Intrusion  Detection  And 
recovery  Scheme  based  on  Assert 

WL/AA  10/10/98  12/31/98  S25000.00  S2396.00 

Image  Quality  Assessment  for  ATR  Applications 
Using  Multiresolurional  Informatio 

WL/ML  01/01/99  12/31/99  $25000.00  $5765.00 

Synthesis  of  New  Two- Photon  Absorbing  Dyes, 

Monomers  and  Polymers 

WL/FI  01/01/99  12/31/99  $25000.00  $7639.00 

Acoustic  Disturbance  Source  Modeling  and 
Development  for  Hypersonic  Receptivity 

WL/ML  01/01/99  12/31*99  $25000.00  $28010.00 

Biocatalysis  of  Biphenyl  and  Diphenylacetylene 
in  an  Aqueous -Organic  Biphasic  Re 

WL/FI  01/01/99  12/3199  $25000.00  $10075.00 

Incorporating  Fixed,  Adaptive,  &  Learning 
Controllers  to  the  Flight  Control 

WL/FI  01/01/99  12/3199  $25000.00  $17448.00 

Dynamics  and  Control  cf  Spacecraft  Formation 
Flying 

WL/FI  01/01/99  12/3199  $25000.00  $10310.00 

Micro-Scale  Flow  Field  Measurement  of  the  Thin 
Meniscus  of  Capillary-Driven  Heat 

WL/AA  01/01/98  12/3198  $25000.00  $13183.00 

Uncertainty  Modeling  cf  Target  Locations  From 
Multiplatform  and  Multisensor  Data 

WL/MN  01/01/99  12/3199  $25000.00  $1991.00 

Sensor  Fusion  w/Passive  Millimeter  Wave  &  Laser 
Radar  for  Target  Detection 

WL/FI  01/01/99  12/3199  $25000.00  $12521.00 

Boundary  Conditions  in  Curvilinear  Coordinates 
for  Direct  Numerical  Simulation 
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SREP  SUB-CONTRACT  DATA 


Report  Author 
Author's  University 


Author's 

Degree 


Sponsoring  _  _  «...  Contract 

Lab  Performance  Period 


Univ.  Cost 
Share 


Mungan  ,  Carl  PhD 

Dept  of  Physics  99-0824 

University  of  Florida,  Pensacola,  FL 


WUMN  01/01/99  12/31/99  524914.00  S3276.00 

infrared  Spectropolarimetr:c  Directonal 
Reflectance  and  Emissivity  zz  Mental  Sur 


Ogale ,  Ainod  PhD 

Chemical  Engineering  99-0917 

Clem  son  University,  Clemson,  SC 


WL/ML  01/01/99  12/31/99  S25000.00  $9000.00 

Structural  Changes  in  Mesophase  Pitch-Based 
Carbon  Fibers*. In  SITU  &ES  SITU  Measu 


Pidaparti ,  Ram  an  a  PhD 

Aeronautics  &  Astronautics  99-0822 

Indiana  U-Purdue  at  Indianap,  Indianapolis,  IN 


WL/FI  01/01/99  12/31/99  $25000.00  $1058100 

Benchmarking  Aerodynamic  Panel  Methods  for 
Flight  Loads  in  Multidisciplinary  Opt 


Saddow ,  Stephen  PhD 

Electrical  Engineering  99-C827 

Mississippi  State  University,  Mississippi  State, 


WIVPO  01/01/98  12/31/98  $25000.00  $0.00 

Silicon  Carbide  Implant  Activation  &  Surface 
preparation  Investigation 


Sepri  ,  Paavo  PhD 

Engineering  Science  99-06  36 

Florida  Inst  of  Technology  ,  Melbourne,  FL 


WL/PO  01/01/99  12/31/99  $25000.00  $6519.00 

Computational  Study  of  Unsteady  Flow 
Interactions  Between  Turbine  31ades,  Cylind 


Shi ,  Hongchi 

Computer  Engineering  99-0825 

Univ  of  Missouri  -  Columbia,  Columbia,  MO 


WL/MN  01/01/99  12/31/99  $25000.00  $15851.00 

Developing  an  efficient  Algorithm  for  Routing 
Processors  of  the  VGI  Parallel  Com 


Soumekh  ,  Mehrdad 
Elec/Computer  Engineering 
SUNT  Buffalo,  Amherst,  NY 


PhD  WL/AA  01/01/99  12/30/99  $25000.00  $0.00 

99-0807  Signal  and  Image  Processing  for  FOPEN/GPEN  SAR 


Rivielk) ,  Craig  BS 

Mechanical  99-0837 

Wright  State  University,  Dayton,  OH 


WL/ML  01/01/99  01/01/99  525000.00  $6304.00 

In-Situ  Synthesis  of  Dis continuously  Reinforced 
Titanium  alloy  Composites  Via  B1 
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APPENDIX  1: 


SAMPLE  SREP  SUBCONTRACT 
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AIR  FORCE  OFFICE  OF  SCIENTIFIC  RESEARCH 
1999  SUMMER  RESEARCH  EXTENSION  PROGRAM 
SUBCONTRACT  99-0806 

BETWEEN 


Research  &  Development  Laboratories 
5800  Uplander  Way 
Culver  City,CA  90230-6608 

.AND 


Penn  State  Erie,  The  Behrend  College 
Contracts  and  Grants  Office 
Erie,  PA  16563 


REFERENCE:  Summer  Research  Extension  Program  Proposal  98-0029 

Start  Date:  01/01/99  End  Date:  12/31/99 

Proposal  Amount:  $25,000 

Proposal  Title:  Image  Quality  Assessment  for  ATR  Applications  Using 

Multiresolutional  Informational  Information  Metrics 


PRINCIPAL  INVESTIGATOR:  Dr.  Farid  Ahmed 

Penn  State  Erie,  The  Behrend  College 
Erie,  PA  16563 

(2)  UNITED  STATES  AFOSR  CONTRACT  NUMBER:  F49620-93-C-0063 

(3)  CATALOG  OF  FEDERAL  DOMESTIC  ASSISTANCE  NUMBER  (CFDA):  12.800 
PROJECT  TITLE:  AIR  FORCE  DEFENSE  RESEARCH  SOURCES  PROGRAM 

(4)  ATTACHMENTS 

1  REPORT  OF  INVENTIONS  .AND  SUBCONTRACT 

2  CONTRACT  CLAUSES 

3  FINAL  REPORT  INSTRUCTIONS 

***SIGN  SREP  SUBCONTRACT  AND  RETURN  TO  RDL*** 
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AIR  FORCE  OFFICE  OF  SCIENTIFIC  RESEARCH 
1999  SUMMER  RESEARCH  EXTENSION  PROGRAM 
SUBCONTRACT  99-0806 
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Research  &  Development  Laboratories 
5800  Uplander  Way 
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Penn  State  Erie,  The  Behrend  College 
Contracts  and  Grants  Office 
Erie,  PA  16563 


REFERENCE:  Summer  Research  Extension  Program  Proposal 

Start  Date:  01/01/99  End  Date: 

Proposal  Amount:  $25,000 

Proposal  Title:  Image  Quality  Assessment  for  ATR  Applications  Using 

Multiresolutional  Informational  Information  Metrics 


PRINCIPAL  INVESTIGATOR:  Dr.  Farid  Ahmed 

Penn  State  Erie,  The  Behrend  College 
Erie,  PA  16563 

(2)  UNITED  STATES  AFOSR  CONTRACT  NUMBER:  F49620-93-C-0063 

(3)  CATALOG  OF  FEDERAL  DOMESTIC  ASSISTANCE  NUMBER  (CFDA):  12.800 
PROJECT  TITLE:  Am  FORCE  DEFENSE  RESEARCH  SOURCES  PROGRAM 

(4)  ATTACHMENTS 

1  REPORT  OF  INVENTIONS  AND  SUBCONTRACT 

2  CONTRACT  CLAUSES 

3  FINAL  REPORT  INSTRUCTIONS 

***SIGN  SREP  SUBCONTRACT  AND  RETURN  TO  RDL*** 


98-0029 

12/31/99 
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1 .  BACKGROUND:  Research  &  Development  Laboratories  (RDL)  is  under  contract 


(749620-93-C-0063)  to  the  United  States  Air  Force  to  administer  the  Summer 
Research  Program  (SRP),  sponsored  by  the  Air  Force  Office  of  Scientific  Research 
(AFOSR),  Bolling  Air  Force  Base,  D  C.  Under  the  SRP,  a  selected  number  of  college 
faculty  members  and  graduate  students  spend  part  of  the  summer  conducting  research 
in  Air  Force  laboratories.  .After  completion  of  the  summer  tour  participants  may 
submit,  through  their  home  institutions,  proposals  for  follow-on  research.  The  follow- 
on  research  is  known  as  the  Summer  Research  Extension  Program  (SREP). 
Approximately  61  SREP  proposals  annually  will  be  selected  by  the  Air  Force  for 
funding  of  up  to  $25,000;  shared  funding  by  the  academic  institution  is  encouraged. 
SREP  efforts  selected  for  funding  are  administered  by  RDL  through  subcontracts  with 
the  institutions.  This  subcontract  represents  an  agreement  between  RDL  and  the 
institution  herein  designated  in  Section  5  below. 

2  RDL  PAYMENTS:  RDL  will  provide  the  following  payments  to  SREP  institutions: 

•  80  percent  of  the  negotiated  SREP  dollar  amount  at  the  start  of  the  SREP 
research  period. 

•  The  remainder  of  the  funds  within  30  days  after  receipt  at  RDL  of  the 
acceptable  written  final  report  for  the  SREP  research. 

3.  INSTITUTION’S  RESPONSIBILITIES:  As  a  subcontractor  to  RDL,  the  institution 
designated  on  the  title  page  will: 
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a.  .Assure  that  the  research  performed  and  the  resources  utilized  adhere  to  those 
defined  in  the  SREP  proposal. 

b.  Provide  the  level  and  amounts  of  institutional  support  specified  in  the  SREP 
proposal.. 

c.  Notify  RDL  as  soon  as  possible,  but  not  later  than  30  days,  of  any  changes  in 
3a  or  3b  above,  or  any  change  to  the  assignment  or  amount  of  participation  of 
the  Principal  Investigator  designated  on  the  title  page. 

d.  .Assure  that  the  research  is  completed  and  the  final  report  is  delivered  to  RDL 
not  later  than  twelve  months  from  the  effective  date  of  this  subcontract,  but  no 
later  than  December  31,  1998.  The  effective  date  of  the  subcontract  is  one 
week  after  the  date  that  the  institution’s  contracting  representative  signs  this 
subcontract,  but  no  later  than  January  15,  1998. 

e.  Assure  that  the  final  report  is  submitted  in  accordance  with  Attachment  3 . 

f  Agree  that  any  release  of  information  relating  to  this  subcontract  (news 
releases,  articles,  manuscripts,  brochures,  advertisements,  still  and  motion 
pictures,  speeches,  trade  associations  meetings,  symposia,  etc.)  will  include  a 
statement  that  the  project  or  effort  depicted  was  or  is  sponsored  by:  Air  Force 
Office  of  Scientific  Research,  Bolling  AFB,  D  C. 

g.  Notify  RDL  of  inventions  or  patents  claimed  as  the  result  of  this  research  as 

specified  in  Attachment  1 . 

h.  RDL  is  required  by  the  prime  contract  to  flow  down  patent  rights  and  technical 
data  requirements  to  this  subcontract.  Attachment  2  to  this  subcontract 
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contains  a  list  of  contract  clauses  incorporated  by  reference  in  the  prime 
contract. 

4.  All  notices  to  RDL  shall  be  addressed  to: 

RDL  AFOSR  Program  Office 
5800  Uplander  Way 
Culver  City,  CA  90230-6609 

5.  By  their  signatures  below,  the  parties  agree  to  provisions  of  this  subcontract. 


Abe  Sopher  Signature  of  Institution  Contracting  Official 

RDL  Contracts  Manager 

Typed/Printed  Name 


Date  Title 


Institution 


Date/Phone 
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ATTACHMENT  2 
CONTRACT  CLAUSES 


This  contract  incorporates  by  reference  the  following  clauses  of  the  Federal  Acquisition 
Regulations  (FAR),  with  the  same  force  and  effect  as  if  they  were  given  in  full  text.  Upon 
request,  the  Contracting  Officer  or  RDL  will  make  their  full  text  available  (FAR  52.252-2). 

FAR  CLAUSES  TITLE  AND  DATE 


52.202- 1 

52.203- 3 

52.203- 5 

52.203- 6 
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INFRARED  ACTIVITY  OF  LITHIUM-DOPED  SOLID  PARAHYDROGEN: 
A  DIFFUSION  QUANTUM  MONTE  CARLO  STUDY 


Robert  J.  Hinde 
Assistant  Professor 
Department  of  Chemistry 
The  University  of  Tennessee 
Knoxville,  TN  37996-1600 


Abstract 


We  have  developed  computational  tools  for  computing  the  lineshape  and  intensity  of  the  Q^O)  IR 
spectral  features  induced  in  solid  p-H2  crystals  by  S-state  atomic  dopants,  and  have  applied  these  tools  to 
Li-doped  solid  p-H2.  The  integrated  absorption  coefficient  for  the  H2  Q^O)  band  for  Li-doped  solid  p-H2 
is  computed  to  be  a  =  4.1  x  10“ 13  cm3/s  per  Li  impurity.  Our  predicted  Q1(0)  lineshape  contains  rich 
structure  which  may  shed  light  on  the  microscopic  morphology  of  Li  trapping  sites  in  solid  p-H2. 
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INFRARED  ACTIVITY  OF  LITHIUM-DOPED  SOLID  PARAHYDROGEN: 
A  DIFFUSION  QUANTUM  MONTE  CARLO  STUDY 


Robert  J.  Hinde 


1.  Introduction 

Atom-doped  solid  H2  has  recently  been  identified  as  an  additive  which  could  improve  the  performance 
of  rocket  propulsion  systems  utilizing  a  liquid  oxygen-liquid  hydrogen  (LOX/LH2)  propellant  [4].  Specific 
impulse  (7sp)  is  a  standard  performance  metric  for  rocket  propulsion  systems,  which  measures  the  momen¬ 
tum  transferred  to  the  rocket  per  unit  weight  of  propellant  consumed.  For  the  current  state-of-the-art 
LOX/LH2  propellant  system,  7sp  =  389  s  [24].  However,  if  the  hydrogen  propellant  contains  A1  atoms  at 
a  mole  concentration  of  5%,  the  theoretical  Jsp  of  this  propellant  system  increases  by  9%  to  425  s;  if  the 
hydrogen  propellant  contains  5%  B  atoms,  the  theoretical  7sp  jumps  to  470  s,  for  an  improvement  of  21%  [4]. 
Consequently,  Fajardo  and  coworkers  at  the  Air  Force  Research  Laboratory  (AFRL)  have  begun  attempts 
to  synthesize  gram-scale  quantities  of  atom-doped  solid  H2  with  impurity  concentrations  approaching  5%. 
The  success  of  these  attempts  would  represent  a  milestone  for  the  Air  Force  High  Energy  Density  Matter 
program. 

To  guide  and  evaluate  these  synthetic  efforts,  the  AFRL  cryosolids  group  must  be  able  to  measure 
accurately  the  concentration  of  atomic  dopants  in  solid  H2  matrices.  Some  of  the  first  experimental  [9, 
10]  and  theoretical  [8]  attempts  to  characterize  atom-doped  H2  matrices  focused  on  the  optical  absorption 
spectra  of  the  dopant  atoms.  As  an  example,  Fig.  1  (reprinted  from  Ref.  [9])  shows  the  absorption  spectrum 
of  solid  H2  doped  with  Li  atoms  at  a  concentration  of  roughly  2  x  10“ 5 .  The  strong  absorption  feature  near 
A  =  670  nm  is  associated  with  the  2s  — Y  2p  electronic  transition  of  the  Li  dopant  atoms. 

Although  electronic  spectroscopy  of  doped  H2  matrices  provides  a  convenient  method  for  detecting  the 
presence  or  absence  of  dopants,  it  is  less  suitable  for  measuring  dopant  concentrations  because  atom-doped 
solid  H2  becomes  optically  thick  in  the  visible  and  ultraviolet  regions  at  very  low  dopant  concentrations. 
For  example,  Fig.  1  indicates  that  the  transmittance  ^670  at  670  nm  for  this  Li-doped  solid  H2  sample  is 
0.8.  A  simple  analysis  based  on  Beer’s  Law  suggests  that  ^670  for  a  sample  of  the  same  thickness  but  with 
a  Li  concentration  of  0.05%  will  be  (0.8)25  «  0.04.  Hence  samples  with  Li  concentrations  of  1  to  5%  will 
be  optically  thick  at  670  nm,  so  that  no  quantitative  concentration  information  can  be  obtained  from  the 
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absorption  spectrum. 

Atom-doped  solid  H2  matrices  also  absorb  weakly  in  the  infrared  (IR)  spectral  region  [11],  as  Fig.  2 
illustrates.  As  we  explain  in  more  detail  below,  the  IR  absorption  features  shown  in  this  figure  arise  from 
nominally  IR  inactive  vibrational  transitions  of  individual  H2  molecules  in  the  matrix  that  become  IR  active 
in  the  presence  of  dopants.  These  features  are  rather  weak,  so  that  1-mm  thick  samples  of  doped  solid  H2  will 
remain  optically  thin  (with  transmittance  values  of  0.1  or  higher)  even  at  dopant  concentrations  approaching 
5%.  Furthermore,  Fig.  2  indicates  that  atom-doped  solid  H2  and  molecule-doped  solid  H2  have  distinct  IR 
spectral  signatures;  this  could  make  it  possible  to  monitor  the  aggregation  of  atomic  dopants  in  H2  matrices 
using  spectroscopic  techniques. 

The  spectra  reprinted  in  Fig.  2  are  some  of  the  first  published  examples  of  infrared  activity  induced  in 
solid  H2  matrices  by  atomic  dopants.  The  integrated  absorption  coefficients  for  the  atom-induced  features 
shown  in  this  figure  are  not  known  experimentally  and  have  not  been  calculated  theoretically.  If  these 
absorption  coefficients  were  known,  spectra  such  as  those  shown  in  Fig.  2  could  obviously  provide  quantitative 
measures  of  the  concentration  of  atomic  impurities  in  doped  solid  H2.  Hence  a  better  understanding  of  these 
spectra  would  be  of  immediate  and  direct  use  to  the  AFRL  cryosolids  group. 

We  report  here  a  calculation  of  the  lineshape  and  absolute  intensity  for  the  IR  spectral  features  induced 
in  solid  H2  by  the  atomic  dopant  Li.  Li-doped  solid  H2  has  been  widely  used  as  a  theoretical  and  experimental 
model  system  in  the  Air  Force  HEDM  program.  The  intensities  of  the  Li-induced  IR  features  are  related 
to  the  total  space-fixed  transition  dipole  moment  vector  induced  in  the  solid  H2  matrix  upon  vibrational 
excitation  of  a  single  H2  molecule  in  the  vicinity  of  the  Li  atom.  We  have  computed  this  transition  moment 
as  the  expectation  value  of  the  corresponding  transition  dipole  operator  over  the  ground  state  wavefunction 
of  the  doped  solid  H2  matrix,  using  diffusion  quantum  Monte  Carlo  (DQMC)  simulations  to  sample  the 
wavefunction  of  the  H2  matrix  and  descendant  weighting  techniques  to  compute  expectation  values  of  the 
transition  dipole  operator. 

The  remainder  of  this  report  is  organized  as  follows.  Section  2  describes  the  theoretical  basis  for 
understanding  the  IR  activity  induced  in  solid  H2  matrices  by  S-state  atomic  dopants.  In  Section  3,  we 
show  how  DQMC  simulations  of  solid  H2  can  be  used  to  compute  the  lineshape  and  intensity  of  these 
dopant-induced  spectral  features.  Section  4  summarizes  our  work  on  Li-doped  solid  H2. 
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2.  Theoretical  background:  IR  activity  induced  in  solid  H2  by  S-state  atomic  dopants 

Here  we  outline  the  theoretical  basis  for  understanding  the  IR  activity  of  solid  H2  matrices  doped  with 
S-state  atomic  impurities.  We  begin  by  summarizing  some  of  the  properties  of  pure  solid  H2,  which  has  been 
the  subject  of  a  number  of  recent  reviews  [23,  29]. 

A.  Properties  of  solid  hydrogen 

In  solid  H2 ,  the  rotational  and  vibrational  motions  of  individual  H2  molecules  are  essentially  unperturbed 
by  the  H2-H2  intermolecular  interactions;  hence  the  rotational  and  vibrational  quantum  numbers  j  and  v  of 
each  molecule  remain  good  quantum  numbers  in  solid  H2  [23,  29].  The  requirements  of  quantum  statistics 
dictate  that  H2  molecules  with  a  total  nuclear  spin  quantum  number  I  =  0  (or  a  singlet  nuclear  spin  state) 
must  have  even  j  values,  while  H2  molecules  with  I  =  1  (or  a  triplet  nuclear  spin  state)  must  have  odd  j 
values.  Molecules  with  7  =  0  (and  even  j )  are  termed  parahydrogen  or  p-H2;  molecules  with  7=1  (and  odd 
j)  are  termed  orthohydrogen  or  o-H2.  Solid  H2  matrices  consisting  of  virtually  100%  p-H2  molecules  with 
7  =  0  can  be  prepared  by  passing  a  flow  of  H2  gas  over  an  ortho  -*  para  conversion  catalyst  and  condensing 
the  resulting  p-H2  molecules  [11]. 

Solid  p-H2  condenses  into  a  hexagonal  close  packed  (hep)  crystal  structure  with  a  nearest-neighbor 
distance  of  7.16  a0  [23].  The  rotational  degrees  of  freedom  of  each  H2  molecule  in  this  crystal  are  described 
by  the  isotropic  spherical  harmonic  Y00(6}<j> );  the  p-H2  crystal  is  thus  composed  of  spherically  symmetric 
objects  with  multipole  moments  that  are  identically  zero,  and  is  held  together  exclusively  by  dispersion  forces. 
The  H2  molecules  in  the  p-H2  crystal  can  be  treated  as  point  particles,  each  with  a  single  vibrational  degree 
of  freedom.  The  translational  zero  point  motion  of  individual  molecules  in  the  solid  H2  matrix  is  rather 
large,  so  that  the  width  of  the  first  peak  in  the  radial  distribution  function  for  solid  H2  is  approximately 
18%  of  the  nearest-neighbor  distance,  even  in  the  T  =  0  K  limit  [23].  Self-diffusion  in  solid  p-H2  is  slow; 
at  T  =  10  K,  a  given  p-H2  molecule  spends  38  ps  at  a  given  lattice  site  before  “hopping”  to  a  neighboring 
site  [23].  Consequently,  it  is  a  good  approximation  to  treat  each  p-H2  molecule  as  permanently  associated 
with  a  specific  crystal  lattice  site  [8]. 

B.  Spectroscopy  of  pure  solid  parahydrogen 

Consider  the  Q^O)  vibrational  excitation  ( v,j )  =  (0,  0)  -4-  (1, 0)  of  a  single  H2  molecule  in  solid  p-H2.* 
This  transition  will  be  IR  active  only  if  the  initial  and  final  states  of  the  crystal  are  connected  by  a  nonzero 

*  In  an  isolated  H2  molecule,  these  two  rovibrational  states  are  separated  by  4161.2  cm"1  of  energy;  in 
solid  H2,  the  Qx(0)  transition  is  slightly  red-shifted  to  4153.0  cm"1  by  H2-H2  intermolecular  interactions  [29]. 
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transition  dipole  moment  matrix  element.  Because  each  H2  molecule  in  the  system  is  spherically  symmetric 
in  both  the  initial  and  the  final  state  of  the  crystal,  there  can  be  no  contribution  to  the  transition  moment 
from  dipoles  induced  on  individual  H2  molecules  by  the  permanent  electrostatic  multipoles  of  neighboring 
molecules.  Hence  the  transition  moment  for  this  vibrational  excitation  arises  from  “overlap  induction”  or 
“exchange  induction”  effects  due  to  the  interpenetration  of  neighboring  molecules’  electron  distributions, 
and  from  the  mutual  polarization  of  neighboring  molecules  due  to  dispersion  interactions. 

The  total  dipole  moment  of  a  pair  of  H2  molecules  with  quantum  numbers  (vj)  =  (0,  0)  must  be  zero 
at  any  H2-H2  distance:  each  individual  H2  molecule  is  spherically  symmetric  by  virtue  of  its  j  =  0  quantum 
number,  and  the  (H2)2  dimer  therefore  has  inversion  symmetry  about  the  dimer’s  center  of  mass.  However, 
if  one  of  the  two  H2  molecules  is  in  the  (u,j)  =  (1,0)  state,  this  inversion  symmetry  is  broken  and  the 
pair  of  molecules  will  possess  a  nonzero  dipole  moment  arising  from  the  exchange  and  dispersion  induction 
mechanisms  outlined  above.  This  makes  the  H2  vibrational  fundamental  IR  active  in  (H2)2  van  der  Waals 
dimers  [16,  17]  and  in  the  collision-induced  absorption  spectrum  of  dense  hydrogen  gas  [19]. 

In  a  pure  p-H2  crystal,  however,  the  QJO)  transition  is  IR  inactive.  Suppose  that  a  particular  H2 
molecule  in  the  hep  p-H2  crystal  is  excited  to  the  (v,j)  =  (1,0)  state.  In  the  hep  lattice,  the  vectors 
connecting  this  excited  H2  molecule  with  its  twelve  nearest  neighbors  sum  to  zero.  Hence  the  net  dipole 
moment  induced  in  the  crystal  by  the  vibrational  excitation  of  one  molecule  is  zero  [28],  although  the 
vibrationally  excited  H2  molecule  induces  “pairwise”  dipole  moments  with  each  of  its  v  =  0  neighbors  as 
shown  in  Fig.  3.  This  consequence  of  the  hep  lattice  symmetry  is  known  as  the  “cancellation  effect”  [27,  28]; 
it  causes  the  Q^O)  transition  to  be  IR  inactive  in  the  pure  p-H2  hep  crystal,  as  is  evident  from  the  bottom 
spectrum  in  Fig.  2. 

C.  Spectroscopy  of  solid  parahydrogen  doped  with  orthohydrogen  molecules 

Now  consider  a  crystal  of  solid  p-H2  in  which  a  single  p-H2  molecule  has  been  replaced  by  an  o-H2 
molecule  with  j  =  1.  The  o-H2  molecule  is  not  spherically  symmetric  and  thus  has  a  nonzero  quadrupole 
moment  (along  with  higher  multipole  moments).  The  quadrupole  field  of  the  o-H2  molecule  induces  dipole 
moments  in  nearby  p-H2  molecules;  the  Qx(0)  transitions  of  these  p-H2  molecules  become  IR  active  by  virtue 
of  these  induced  moments. 

To  demonstrate  this,  we  compute  the  transition  dipole  moment  associated  with  this  excitation  using 
arguments  presented  in  Ref.  [21],  making  the  simplifying  assumption  that  the  centers  of  mass  of  the  H2 
molecules  in  the  crystal  are  held  fixed.  Let  En  be  the  electric  field  of  the  o-H2  impurity  at  the  nth  p-H2 
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molecule;  the  dipole  moment  induced  in  this  p-H2  molecule  by  the  o-H2  impurity  is  mn  =  47re0aEn  where 
a  is  the  (isotropic)  polarizability  of  the  p-H2  molecule.  (We  neglect  the  overlap  induced  dipole  described 
above;  at  the  nearest-neighbor  distance  of  solid  H2,  this  dipole  is  much  smaller  in  magnitude  than  the  dipole 
induced  by  the  quadrupole  field.)  If  we  assume  that  the  total  induced  dipole  M  of  the  system  is  simply  the 
sum  of  these  induced  moments,  then  the  transition  moment  associated  with  the  Q^O)  transition  of  the  Arth 
p-H2  molecule  is 

=  0  4™oa  Yj  En  f)  =  4jrco°oiE, fc  (1) 

where  |i)  and  |/)  represent  the  initial  and  final  states  of  the  doped  crystal,  the  superscript  ( k )  indicates 
that  the  kth.  p-H2  molecule  is  undergoing  excitation,  and  a01  =  (y  =  0|or|v  =  l)  is  the  isotropic  “transition 
polarizability”  of  H2,  which  has  the  value  or01  =  0.65  au  [20].  Note  that  terms  in  M  with  n  ^  k  do  not 
contribute  to  due  to  the  orthogonality  of  the  H2  vibrational  wavefunctions  with  vk  =  0  and  vk  =  1. 
Because  ^  0  if  ^  0,  the  Q:(0)  transition  of  a  p-H2  molecule  near  an  o-H2  impurity  will  be  IR 

active.  The  top  spectrum  shown  in  Fig.  2  is  the  IR  absorption  spectrum  of  a  p-H2  matrix  doped  with  2% 
o-H2;  the  Qi(0)  transition  of  the  p-H2  host  is  clearly  visible. 

The  integrated  absorption  coefficient  of  this  transition  is  related  to  through  [1] 


Although  this  sum  appears  to  diverge  as  the  number  of  p-H2  molecules  JV  ->  oo,  the  contribution  to  a 
from  a  p-H2  molecule  which  is  far  away  from  the  o-H2  impurity  is  negligible  because  E^  decays  quickly  with 
increasing  distance.  Evaluation  of  using  Eqs.  (1)  and  (2)  leads  to  a  theoretical  Qx(0)  integrated  IR 
absorption  coefficient  ofor  =  2.0xl0”14  cm3/s  per  o-H2  impurity  [21];  the  corresponding  experimental  value 
is  2.2  x  10” 14  cm3/s  [12].  The  small  discrepancy  between  theory  and  experiment  can  be  rationalized  by 
recalling  that  we  have  neglected  both  the  motion  of  the  H2  molecules  and  the  contributions  to  the  transition 
dipole  moment  from  overlap  induction  and  from  three-body  effects. 

A  similar  approach  has  recently  been  used  to  estimate  the  integrated  absorption  coefficients  of  the  S1(0) 
and  U1(0)  transitions  in  solid  p-H2  doped  with  o-H2  impurities  [2];  the  good  agreement  between  experiment 
and  theory  obtained  for  these  transitions  underscores  the  validity  of  this  approach. 


D.  Spectroscopy  of  solid  parahydrogen  doped  with  S-state  atoms 

Finally  we  consider  a  p-H2  crystal  with  N  H2  molecules  and  a  single  S-state  atomic  impurity.  As  Fig.  2 
shows,  S-state  atoms  embedded  in  p-H2  induce  weak  IR  activity  in  the  Qx(0)  spectral  region.  In  contrast 
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to  o-H2  impurities,  S-state  atoms  have  no  permanent  multipole  moments;  hence  the  transition  moment 
associated  with  this  IR  activity  must  arise  solely  from  overlap  and  dispersion  effects,  and  from  the  fact  that 
the  dopant  atom  breaks  the  symmetry  of  the  hep  lattice.  We  now  derive  an  expression  for  this  transition 
dipole  moment. 

Let  R0  denote  the  position  of  the  atomic  impurity  and  {Rn:n=l,2,3,...,7V}  denote  the  centers  of 
mass  of  the  H2  molecules  in  the  crystal.  The  vibrational  quantum  number  of  the  H2  molecule  at  Hk  will  be 
denoted  vk ;  we  assume  that  the  rotational  quantum  number  j  of  each  H2  molecule  is  zero.  (The  quantum 
number  j  will  remain  a  good  quantum  number  if  the  dopant~H2  interaction  is  isotropic,  or  very  nearly  so,  at 
low  energies;  this  is  the  case  for  Li  dopants  [6].)  The  state  of  the  doped  crystal  is  therefore  completely  specified 
by  the  position  vectors  {  R0 ,  Kx ,  R2 , . . . ,  R^  }  and  the  list  of  vibrational  quantum  numbers  ,  ^2’  •  -  •  >  }  * 

We  focus  our  attention  on  the  vx  =  0  ->  1  transition  of  the  H2  molecule  at  R:.  Let  | i)  denote  the  initial 
state  of  the  doped  crystal,  in  which  vk  =  0  for  all  k ,  and  let  |/)  denote  the  final  state  of  the  crystal,  in  which 
vx  =  1  and  vk  =  0  for  k  >  1.  We  assume  that  the  dipole  moment  operator  for  this  system  can  be  written  as 

a  pairwise  sum  of  overlap-induced  dipoles: 

N- 1  N  N  N 

M(1)  =  mU,k)  = 

j= o  k=j+ 1  k~l  k  = 2 


Em(0'fc)+E 


m 


(M) 


(3) 


where  represents  the  dipole  moment  induced  by  the  overlap  of  the  particles  at  R j  and  R^.  Note  that 

if  the  dipole  moment  is  pairwise  additive,  (p-H2)2  pairs  in  which  both  molecules  remain  in  the  ground  state 
(v}j)  =  (0,0)  cannot  contribute  to  the  dipole  moment. 

The  transition  moment  for  the  Q:(0)  transition  of  the  H2  molecule  at  Rx  is  =  <i|ivr(1)|/).  To 
evaluate  ,  we  first  integrate  over  the  vibrational  coordinates  of  all  N  H2  molecules.  Because  the  v  =  0 
and  v  =  1  vibrational  wavefunctions  for  H2  are  orthogonal,  this  integration  eliminates  from  MyV  all  terms 


which  do  not  depend  on  the  vibrational  coordinate  of  the  molecule  at  Rx.  Hence  we  obtain 

N 


mJV  =  J  »J(Q) 


m 


(o,i) 

01 


+E 


m 


(i^) 

01 


k  —  2 


^y(Q)dQ 


(4) 


where  Q  =  {R0,  R1?  R2, . . . ,  HN}  is  the  collection  of  position  vectors  describing  the  system.  The  initial  and 
final  spatial  wavefunctions  of  the  system  (after  integrating  over  the  H2  vibrational  coordinates)  are  ^-(Q) 
and  ^(Q).  The  vibrationally  averaged  pairwise  contributions  to  the  transition  moment  are 

m£ll},1)  =  (vx  =  Olmt0-1^!  =  l)  (5) 


for  the  dipole  arising  from  overlap  between  the  atomic  impurity  and  the  H2  molecule  at  R1?  and  (for  k  >  1) 


=  (t>!  =  0,vk  =  =  l,vk  =  0) 


(6) 
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for  the  dipole  arising  from  overlap  between  the  H2  molecules  at  Rx  and  Kk.  We  emphasize  that  mg^  and 
are  vector  quantities  which  are  oriented  along  the  vectors  R01  =  Rx  —  R0  and  Rx  k  =  Hk  —  R: 
respectively,  and  which  depend  on  the  distances  Ro  x  and  Rl  k . 

From  Eq.  (4)  we  obtain  the  total  transition  dipole  moment  for  vibrational  excitation  of  the  H2  molecule 
located  at  Rx.  The  contribution  to  the  crystal’s  integrated  absorption  coefficient  from  this  excitation  is 
[cf.  Eq.  (2)] 

=  (7) 

Each  H2  molecule  contributes  in  a  similar  fashion  to  the  absorption  spectrum,  so  that  the  total  absorption 


coefficient  per  atomic  impurity  is 


4-?*B,-S3EX>r 


(8) 


3.  Computational  approach:  DQMC  simulations  of  the  H2  Q1(0)  band 


A.  Absorption  coefficients 

The  transition  moment  is  an  integral  over  the  initial  wavefunction  ^(Q)  and  the  final  wavefunction 
\£y(Q)  of  the  doped  p-H2  matrix.  Because  solid  p-H2  is  a  quantum  crystal,  individual  H2  molecules  in  the 
crystal  undergo  large-amplitude  anharmonic  zero  point  motion  around  their  nominal  lattice  sites  [23].  In 
addition,  the  instantaneous  positions  of  neighboring  p-H2  molecules  are  highly  correlated  [8].  Therefore  the 
many-body  wavefunctions  and  Wy  characterizing  the  doped  crystal’s  translational  degrees  of  freedom 
cannot  be  obtained  in  analytic  form,  and  Eq.  (4)  cannot  be  evaluated  directly.  However,  we  can  estimate 
using  Monte  Carlo  sampling  techniques. 

Without  loss  of  generality,  we  can  take  ^(Q)  to  be  a  real  function,  so  that  ^,*(Q)  =  ^(Q).  Then 


Eq.  (4)  can  be  rewritten  as 


M$V  =  /*i(Q)  +  ^(Q)dQ  =  J  /(Q)P(Q)dQ 


where 


/(  q) =mr+£x 


is  the  instantaneous  Qx(0)  transition  moment  for  crystal  configuration  Q,  and  P(Q)  =  ^(Q)^  (Q)  can  be 
interpreted  as  a  probability  density.  If  we  have  a  method  for  selecting  configurations  {Qi,Q2jQ3>  —  *}  from 
P( Q),  then  we  can  approximate  the  integral  as 


mS?«i/»x;/(q,). 

j=i 
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As  the  number  of  sampled  configurations  n  oo,  this  Monte  Carlo  approximation  converges  to  the  integral’s 
exact  value. 

There  is  no  computationally  efficient  method  for  selecting  configurations  according  to  a  probability 
density  of  the  form  ${(Q)$y(Q)  when  the  wavefunctions  and  \ty  are  both  unknown.  However,  if  we 
make  the  simplifying  assumption  that  \Pf.  «  \Py ,  then  we  can  replace  P( Q)  in  Eq.  (9)  with  P(Q)  =  |^t-(Q)|2; 
a  probability  density  function  with  this  form  can  be  readily  sampled  during  a  DQMC  simulation*  of  the 
doped  p-H2  crystal  using  descendant  weighted  DQMC  techniques  [5].  This  permits  us  to  estimate  My*\ 
In  fact,  using  descendant  weighting  techniques,  we  can  simultaneously  evaluate  My^  for  each  H2  molecule 
1  <  j  <  N  in  one  DQMC  simulation,  thereby  obtaining  all  of  the  information  needed  to  compute  a  using 
Eq.  (8). 

This  approach  rests  on  the  assumption  that  «  tfy.  The  validity  of  this  assumption  would  be 
questionable  if  vibrationally  excited  p-H2  could  not  easily  fit  into  the  solid  H2  lattice;  in  this  case,  the 
vibrationally  excited  molecule  would  experience  substantial  repulsive  forces,  and  we  would  expect  the  Q^O) 
transition  of  this  molecule  to  be  strongly  blue-shifted.  However,  as  Fig.  2  shows,  the  H2  Qx(0)  transition  in 
doped  solid  H2  is  weakly  red-shifted  from  the  gas-phase  Q^O)  transition  of  4161.2  cm-1;  this  suggests  that 
the  assumption  that  «  \£y  is  reasonable.** 

The  validity  of  our  approach  is  demonstrated  by  our  recent  calculation  [13]  of  the  absorption  coefficient 
for  the  Qx(0)  +  Q2(0)  double  vibrational  transition  in  solid  p-H2  observed  by  Mengel  and  coworkers  [18]. 
In  this  transition,  a  single  IR  photon  excites  the  vibrational  coordinates  of  two  neighboring  p-H2  molecules 
simultaneously;  the  intensity  of  the  transition  arises  from  overlap-induced  (p-H2)2  dipole  moments  described 
by  Eq.  (6).  Because  these  moments  vary  strongly  with  distance  [19],  we  must  explicitly  incorporate  the 
molecules’  translational  zero-point  motion  into  our  calculation  of  the  overall  transition  moment  My,-. 

Using  the  techniques  described  in  this  section,  we  obtained  [13]  an  absorption  coefficient  for  the 

*  DQMC  simulations  sample  the  ground  state  wavefunction  of  the  doped  p-H2  crystal,  and  thus  are 
strictly  valid  only  in  the  limit  T  =  OK.  The  methods  outlined  in  this  proposal  could  also  be  applied  to 
finite-temperature  path  integral  Monte  Carlo  studies  of  solid  p-H2.  However,  thermal  effects  on  the  solid 
p-H2  induced  Qx(0)  intensity  and  lineshape  are  negligible  below  T  =  10  K  [15],  so  a  simpler  DQMC-based 
approach  is  sufficient. 

**  We  note  that  this  assumption  could  be  tested  by  computing  the  transition  moment  My^  using  “nested 
walk”  DQMC  techniques  originally  designed  for  the  quantum  Monte  Carlo  evaluation  of  electronic  transition 
dipole  moments  [3].  However,  these  techniques  are  very  computationally  intensive. 
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Qi(0)  +  Q2(0)  double  transition  of  a  =  1.1  x  10"18  cm3/s;  the  experimental  value  [18]  is  a  =  (9.3  ± 
0.9)  x  1CT19  cm3/s.  If  we  make  the  simplifying  assumption  that  the  H2  molecules  are  fixed  at  their  nom¬ 
inal  lattice  sites,  we  instead  obtain  a  theoretical  absorption  coefficient  of  a  =  2.5  x  10~2°  cm3/s,  which  is 
nearly  40  times  too  small.  This  discrepancy  emphasizes  that  computational  methods  which  account  for  the 
quantum  behavior  of  the  p-H2  crystal,  such  as  DQMC,  are  needed  for  a  quantitative  treatment  of  the  IR 
transitions  in  solid  p-H2  arising  from  overlap-induced  dipole  moments. 

B.  Lineshapes 

The  integrated  absorption  coefficient  a  in  Eq.  (8)  obtained  from  the  individual  transition  moments 
represents  the  total  absorption  intensity  (per  dopant)  of  all  of  the  dopant-induced  IR  spectral  features 
associated  with  the  H2  Qi(0)  transition.  Figure  2  shows  that  molecular  dopants  induce  multiple  absorption 
features  in  the  H2  Q^O)  spectral  region,  and  that  both  atom-  and  molecule-induced  features  apparently 
have  a  finite  width.  The  integrated  absorption  coefficient  a  of  Eq.  (8)  thus  represents  an  integral  over  the 
entire  induced  Q^O)  band,  which  extends  across  the  interval  4140  cm”1  to  4170  cm-1. 

The  multiplicity  of  features  induced  by  molecular  dopants  is  presumably  related  to  the  anisotropy 
of  these  dopants;  p-H2  molecules  on  different  “sides”  of  the  dopants  may  absorb  IR  radiation  of  slightly 
different  frequencies.  To  understand  how  dopant-induced  IR  spectra  of  solid  H2  could  be  used  to  monitor 
the  aggregation  of  dopants  in  the  matrix,  we  must  therefore  be  able  to  predict  the  lineshape  of  the  Qi(0) 
band  as  well  as  its  total  intensity.  We  now  outline  a  method  for  predicting  this  lineshape  using  DQMC 
simulations. 

To  perform  a  DQMC  simulation  of  doped  solid  p-H2,  we  require  a  potential  energy  surface  for  the 
doped  p-H2  crystal.  This  surface  is  a  function  of  the  positions  of  the  individual  p-H2  molecules  and  the 
dopant,  and  of  the  vibrational  quantum  numbers  of  the  p-H2  molecules:  V  =  V(Qyvvv2, . . .  ,vN).  A 
DQMC  simulation  of  the  crystal  under  the  constraint  vx  =  v2  =  •  *  •  =  vN  =  0  would  sample  configurations 
from  the  initial  wavefunction  tft-(Q)  and,  if  coupled  with  descendant  weighting  techniques,  would  allow  us 
to  compute  expectation  values  over  |^t-(Q)|2.  The  DQMC  simulation  also  generates  an  estimate  of  the  total 
energy  E0  of  the  ground  state  of  the  doped  crystal  under  the  constraint  =  v2  =  •  •  •  =  =  0. 

If  we  perform  a  second  DQMC  simulation  of  the  crystal,  this  time  under  the  constraint  vx  =  I  and 
v2  =  v3  =  •  •  •  =  vN  =  0,  we  obtain  the  energy  e[^  of  the  state  of  the  crystal  in  which  molecule  1  has  been 
excited  to  the  v  =  1  level.  (We  assume  for  the  moment  that  the  vibrational  excitation  remains  localized 
on  molecule  1.)  The  transition  energy  between  the  ground  state  of  the  doped  crystal  and  this  vibrationally 
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excited  state  is  A =  E[ ^  —  E0 .  According  to  this  picture,  the  Qj(0)  transition  of  this  p-H2  molecule 
would  be  a  sharp  peak  at  this  transition  energy,  with  a  total  intensity  given  by  Eq.  (7). 

However,  if  the  vibrational  coordinate  of  molecule  1  is  only  weakly  coupled  to  the  translational  coor¬ 
dinates  Q,  which  is  to  say  that  ^(Q,vx  =  0,  v2, . . . ,  vN)  «  ^(Q,^  =  1,  v2} . . . ,  vN),  then  AE^  can  be 
approximated  using  perturbation  theory.  The  ground  state  energy  E0  satisfies 

=  <®i  |^|**>  =  <®i  I#"  +  ^(Qf  =  0f  - . . ,  (12) 

and  E[1>}  satisfies 

£<1}  =  <¥,!%,)  =  (^|f  +V(Q,Vl  =  1,«2, (13) 

But  if  we  define 

A^Q)  =  V(Q,  Wl  =  1,  v2 . vN)  -  V(Q,  Vl  =  0,v2,...,  vN),  (14) 

and  if  AV'W(Q)  is  small  at  configurations  Q  which  are  accessible  to  the  ground  state  wavefunction  (which 

implies  «  \Py),  then  from  first-order  perturbation  theory  we  have 

e[1]  «  £*0+  («<|AV’^(Q)|$1-)  or  A E™  »  1 AK(1) (Q) | > .  (15) 

This  expectation  value  of  AV^1)  over  |^,*(Q)|2  can  be  computed  using  descendant  weighting  techniques 

during  the  same  DQMC  simulation  in  which  we  compute  Therefore  with  a  single  DQMC  simulation, 

we  can  compute  the  transition  energies  A E^  and  absorption  coefficients  for  each  p-H2  molecule  in  the 
doped  crystal. 

The  simulated  spectrum  which  would  result  from  such  a  calculation  would  be  a  stick  spectrum,  composed 
of  sharp  peaks  at  the  transition  energies  A E^  with  respective  absorption  coefficients  The  spectra 

shown  in  Fig.  2  are  clearly  not  stick  spectra;  the  dopant-induced  Qx(0)  features  have  finite  widths.  Although 
vibron  hopping  in  the  p-H2  crystal  broadens  IR  spectral  features  in  solid  H2  [29],  recent  high-resolution 
IR  studies  [7,  15]  demonstrate  that  this  effect  is  relatively  weak,  and  broadens  features  in  the  Qx(0)  band 
only  by  about  0.1  cm-1.  Consequently  it  is  a  good  approximation  to  assume  that  the  vibrational  excitation 
remains  localized  on  a  single  p-H2  molecule.  The  widths  of  the  peaks  shown  in  Fig.  2  therefore  seem  to  arise 
primarily  from  instrumental  broadening;  the  resolution  of  the  spectrometer  used  to  record  these  spectra  is 
1  cm-1.  We  emphasize  that  any  broadening  due  to  vibron  hopping  will  not  change  the  integrated  absorption 
coefficient  of  the  Q^O)  band;  this  is  a  consequence  of  the  principle  of  spectroscopic  stability  [27]. 
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4.  Studies  of  Li-doped  solid  H2 


A.  Quantum  chemical  calculations 

The  DQMC  simulations  described  above  require  as  input  both  (1)  the  potential  energy  surface  V(Q,  vv 
v2, . . . ,  vN)  of  the  doped  p-H2  crystal  and  (2)  induced  transition  moment  functions  for  Li-p-H2  and  (p-H2)2 
pairs.  We  assume  that  V  is  pairwise  additive,  and  that  the  interaction  between  two  v  =  0  p-H2  molecules 
is  described  by  the  Silvera-Goldman  [22]  potential  function  VSG(J?);  this  potential  includes,  in  an  effective 
form,  three-body  contributions  to  the  potential  energy  which  are  important  in  solid  H2.  To  describe  the 
interaction  between  a  p-H2  molecule  with  v  =  0  and  a  p-H2  molecule  with  v  =  1,  we  add  to  VSG  ( R)  an 
incremental  contribution  obtained  from  Ref.  [19],  which  is  the  difference  between  the  (vl}v2)  =  (0,0)  and 
{v11v2)  =  (0, 1)  potential  curves  for  an  isolated  gas-phase  (p-H2)2  dimer. 

The  Li-(p-H2)  interaction  was  computed  ab  initio  at  the  counterpoise-corrected  CCSD(T)  level  [26], 
using  aug-cc-pVTZ  basis  sets  on  the  hydrogen  atoms,  a  truncated  cc-pVQZ  basis  set  (which  omits  g  orbitals) 
on  Li,  and  an  uncontracted  (3s3p2d)  set  of  bond  functions  midway  between  Li  and  H2.  (Our  results  did 
not  change  significantly  when  either  diffuse  functions  or  tight  core-valence  functions  were  added  to  the  Li 
basis  set.)  Single  point  ab  initio  calculations  of  the  Li-H2  potential  energy  were  performed  at  369  geometries 
using  standard  Jacobi  coordinates;  these  were  averaged  over  the  appropriate  H2  vibrational  wavefunction 
and  the  isotropic  Y0 o(0,  <t>)  spherical  harmonic  to  obtain  Li-(p-H2)  potential  curves  for  a  p-H2  molecule  with 
the  vibrational  quantum  number  v  =  0  or  v  —  1. 

The  H2-H2  induced  transition  moment  function  defined  in  Eq.  (6)  was  taken  directly  from  Ref.  [19]. 
The  Li-H2  induced  transition  moment  function  was  computed  ab  initio  at  the  CCSD(T)  level  [26],  using 
d-aug-cc-pVQZ  basis  sets  on  the  hydrogen  atoms  and  the  truncated  Li  cc-pVQZ  basis  set  described  above; 
no  bond  functions  were  employed.  Single  point  ab  initio  dipole  moments  were  obtained  at  189  geometries, 
again  using  Jacobi  coordinates;  these  dipole  moments  were  averaged  over  the  H2  vibrational  wavefunctions 
and  the  isotropic  Yoo{0}(j>)  spherical  harmonic  to  obtain  the  Li-(p-H2)  induced  transition  moment  function 
defined  in  Eq.  (5). 

Figure  4  shows  our  computed  Li-(p-H2)  potential  energy  and  transition  moment  curves.  Vibrational 
excitation  of  the  p-H2  molecule  strengthens  the  long-range  Li-(p-H2)  attractive  dispersion  interactions  and 
deepens  the  potential  slightly.  A  Numerov-Cooley  analysis  of  the  Li-(p-H2)  bound  states  shows  that  the 
v  =  0  and  v  =  1  Li-(p-H2)  curves  each  support  one  van  der  Waals  bound  state;  the  energies  of  these  van  der 
Waals  complexes  are  E  —  —3.94  cm’1  for  v  =  0  p-H2  and  E  =  —5.12  cm”1  for  v  =  1  p-H2.  The  first-order 
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perturbation  theory  argument  outlined  previously  predicts  that  the  van  der  Waals  bond  should  be  1.10  cm"1 
stronger  in  the  v  =  1  complex.  This  is  in  good  agreement  with  the  Numerov-Cooley  results,  which  show 
that  vibrational  excitation  of  the  p-H2  molecule  strengthens  the  bond  by  1.18  cm-1. 

B.  DQMC  simulations  of  Li-induced  spectral  features 

We  performed  DQMC  simulations  of  Li-doped  solid  p-H2  using  a  protocol  similar  to  that  employed 
by  Cheng  and  Whaley  [8].  We  began  with  a  DQMC  simulation  “box”  containing  180  p~H2  molecules 
arranged  on  the  hep  lattice.  One  p-H2  molecule  was  removed  from  the  simulation  box  and  replaced  by  a  Li 
atom.  We  then  applied  periodic  boundary  conditions  in  three  dimensions  and  conducted  a  series  of  DQMC 
simulations,  using  the  trial  wavefunction  \Ptr  defined  in  Ref.  [8]  to  reduce  the  statistical  fluctuations  inherent 
in  DQMC  methods  and  increase  the  efficiency  of  our  simulations.  This  trial  function  consists  of  a  product 
of  Gaussians  localizing  each  particle  near  its  nominal  lattice  site  and  two-body  Jastrow  functions  correlating 
the  instantaneous  positions  of  neighboring  particles. 

Our  DQMC  simulations  are  conducted  using  an  imaginary  time  step  of  50  i  au  and  typically  consist  of 
1000  to  2500  time  steps,  beginning  from  a  previously  equilibrated  ensemble  of  250  or  500  configurations. 
During  the  simulations,  we  used  descendant  weighting  techniques  to  evaluate  the  transition  moment 
and  transition  energy  A for  the  twelve  p-H2  molecules  which  are  nearest  neighbors  to  the  Li  impurity.  We 
note  that  descendant  weighting  methods  effectively  eliminate  any  trial  function  bias  in  the  computation  of 
expectation  values  [5],  and  are  therefore  superior  to  “mixed”  estimators  of  observables  which  have  been  widely 
used  in  previous  studies  of  solid  H2  [8].  During  the  first  90%  of  each  DQMC  simulation,  we  accumulated 
expectation  values  via  descendant  weighting. 

From  ten  independent  DQMC  simulations,  we  obtained  an  estimate  for  the  total  Qx  (0)  transition  mo¬ 
ment  of  Myt*  =  (2.25  ±  0.14)  x  10"3  au  which  corresponds  to  a  Q1(0)  integrated  absorption  coefficient  of 
5  =  4.1  x  10~13  cm3/s  per  Li  impurity.  To  place  this  absorption  coefficient  in  context,  the  Qx(0)  integrated 
absorption  coefficient  for  the  transition  induced  in  solid  p-H2  by  the  permanent  quadrupole  moment  of  o-H2 
impurities  (shown  in  the  top  spectrum  of  Fig.  2)  is  d  =  2.2  x  10-14  cm3/s  per  o-H2  impurity.  The  Li-induced 
features  are  thus  roughly  18  times  stronger  than  those  induced  by  o-H2.  (Note  that  the  top  spectrum  of 
Fig.  2  was  taken  at  an  o-H2  concentration  of  2%.) 

Figure  5  shows  our  simulated  absorption  spectrum  for  a  1-mm  thick  sample  of  solid  p-H2  doped  with 
Li  at  a  concentration  of  1000  ppm  [14].  This  spectrum  was  obtained  from  the  theoretical  transition  energies 
A and  absorption  coefficients  5^')  of  the  twelve  p-H2  molecules  surrounding  the  Li  dopant;  the  computed 
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stick  spectrum  was  convolved  with  a  Gaussian  of  width  1  cm-1  to  simulate  instrumental  broadening.  We 
also  show  in  Fig.  5  the  experimental  absorption  spectrum  [11]  for  Na-doped  solid  p-H2,  which  we  have  chosen 
as  a  chemically  similar  system  for  which  an  experimental  lineshape  is  available. 

The  comparison  between  our  simulated  spectrum  for  Li-doped  p-H2  and  the  experimental  results  for  Na- 
doped  p-H2  is  striking.  Apart  from  an  overall  shift,  both  spectra  have  very  similar  lineshapes.  Perhaps  most 
surprising  is  the  fact  that  the  lineshapes  are  clearly  asymmetric,  suggesting  that  either  the  dopants  move  to 
one  side  of  the  trapping  site  or  the  trapping  site  undergoes  a  symmetry-lowering  structural  transformation. 
If  the  Li  and  Na  dopants  remained  symmetrically  distributed  in  the  vacancy  created  by  removing  one  p-H2 
molecule  from  the  hep  lattice,  we  would  expect  to  see  at  most  two  features  in  the  Qx(0)  band:  one  from 
p-H2  molecules  in  the  same  hep  crystal  plane  as  the  dopant,  and  one  from  p-H2  molecules  in  adjacent  crystal 
planes.  The  experimental  spectrum  consists  of  at  least  three  features,  however. 

The  Qi(0)  features  induced  by  Na  dopants  are  red-shifted  from  our  theoretical  Li-induced  features  by 
about  6  cm-1.  The  overall  shift  is  controlled  by  the  change  AV  in  the  potential  surface  induced  upon 
excitation  of  a  single  H2  molecule  [Eq.  (14)],  which  is  dominated  by  the  increase  in  attractive  dopant-(p- 
H2)  dispersion  interactions  which  accompanies  H2  vibrational  excitation.  If  we  naively  assume  that  these 
dispersion  interactions  vary  with  the  dopant’s  polarizability,  we  would  expect  the  Na-induced  features  to  be 
more  highly  red-shifted  than  the  Li-induced  features. 

5.  Summary 

We  have  developed  computational  tools  for  computing  the  lineshape  and  intensity  of  the  Q1(0)  IR 
spectral  features  induced  in  solid  p-H2  crystals  by  S-state  atomic  dopants,  and  have  applied  these  tools  to 
the  model  HEDM  system  of  Li  in  solid  p-H2.  The  integrated  absorption  coefficient  for  the  H2  Q1(0)  band  for 
Li-doped  solid  p-H2  is  computed  to  be  d  =  4.1  x  10“13  cm3/s  per  Li  impurity.  Our  predicted  Q1(0)  lineshape 
(Fig.  5)  contains  rich  structure  which  may  shed  light  on  the  microscopic  morphology  of  Li  trapping  sites 
in  solid  p-H2;  this  lineshape  agrees  very  well  with  the  experimental  Qx(0)  lineshape  [11]  for  the  chemically 
related  system  of  Na-doped  solid  p-H2. 
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7.  Figure  captions: 

Fig.  1.  Visible  transmission  spectrum  of  Li-doped  solid  H2.  The  dashed  line  is  the  spectrum  of  a  sample  of 
Li  in  solid  H2  at  a  concentration  of  approximately  20  ppm  and  at  T  =  3.4  K;  the  dotted-dashed  line  is  the 
spectrum  of  pure  solid  H2  at  T  =  3.4  K.  The  solid  curve  is  the  difference  spectrum.  The  sample  thickness  is 
estimated  to  be  roughly  10  /im.  Reprinted  from  Ref.  [9]. 

Fig.  2.  Infrared  absorption  peaks  induced  in  solid  p-H2  samples  by  various  atomic  and  molecular  dopants. 
The  temperature  of  each  sample  is  T  =  2  K;  the  sample  thicknesses  range  from  1  to  3  mm.  Except  for  o-H2 
dopants  (top  spectrum),  the  estimated  dopant  concentrations  range  from  100  to  1000  ppm.  Reprinted  from 
Ref.  [11]. 

Fig.  3.  Illustration  of  the  cancellation  effect  in  hep  solid  p-H2.  The  open  circles  represent  v  =  0  p-H2 
molecules;  the  filled  circle  represents  a  v  =  1  p-H2  molecule.  The  vibrationally  excited  molecule  induces 
pairwise  dipole  moments  with  its  neighbors  (arrows),  but  these  moments  sum  to  zero  due  to  the  symmetry 
of  the  hep  crystal  lattice. 

Fig.  4.  Ab  initio  Li-(p-H2)  potential  energy  curves  (upper  panel)  and  Q^O)  transition  moment  (lower  panel); 
see  the  text  for  computational  details.  In  the  upper  panel,  the  solid  line  indicates  v  =  0  p-H2  and  the  dashed 
line  indicates  v  =  1  p-H2. 
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Fig.  5.  Theoretical  spectrum  for  the  Li-induced  H2  Qx(0)  band  in  Li-doped  solid  p-H2  (lower  spectrum) 
compared  with  the  experimental  spectrum  for  the  Na-induced  H2  Q^O)  band  in  Na-doped  solid  p-H2  (upper 
spectrum,  taken  from  Ref.  [11]).  The  simulated  Li-induced  spectrum  corresponds  to  a  1-mm  thick  sample  of 
Li-doped  solid  H2  with  a  Li  concentration  of  1000  ppm;  the  absorbance  scale  for  this  spectrum  is  absolute. 
The  experimental  Na-induced  spectrum  is  for  a  sample  about  1-mm  thick  with  an  estimated  Na  concentration 
of  100  to  1000  ppm;  the  absorbance  scale  for  this  spectrum  is  not  the  same  as  for  the  simulated  Li-induced 
spectrum. 
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Abstract 

Blind  point-source  image  restoration  refers  to  the  problem  of  high  resolution  recovery  of  point¬ 
like  sources  which  are  blurred  by  an  unknown  point  spread  function  (psf)  and  corrupted  by 
noise.  Applications  include  astronomical  star  field  localization,  airborne  target  image  localiza¬ 
tion,  magnetoencephalogram  brain  current  imaging,  and  seismic  deconvolution.  Both  single 
and  multiple  frame  observation  cases  are  addressed,  with  an  emphasis  on  adaptive  optics  (AO) 
telescope  data  sets.  This  final  report  describes  two  methods  for  solving  this  problem  which 
were  developed  under  AFOSR  SREP  grant  support.  First,  it  will  be  shown  that  with  suit¬ 
able  constraints,  the  problem  can  be  cast  in  the  language  of  subspace  decomposition  as  used  in 
blind  signal  copy  algorithms  for  digital  wireless  communications.  Assuming  a  separable  psf,  we 
propose  a  deterministic,  non-iterative  ESPRIT-like  solution  to  the  restoration  problem.  This 
algorithm  is  then  extended  to  apply  to  non-separable  psf’s  by  approximating  them  as  a  series 
expansion  of  a  few  separable  components. 

The  second  algorithm  introduced  in  this  report  is  a  Bayesian  method  for  blind  restoration  of 
images  of  sparse,  point-like  objects.  The  proposed  method  uses  maximum  a  posteriori  estimation 
techniques  to  recover  both  the  unknown  object  and  blur.  Markov  random  field  (MRF)  models 
are  used  to  represent  prior  information  about  both  the  sparse,  point-like  structure  of  the  object, 
and  the  smoothed  random  structure  of  the  blur.  As  compared  with  general  purpose  blind 
algorithms,  incorporating  a  sparse  point  source  MRF  model  enables  much  higher  resolution 
restorations,  improves  point  localization,  and  aids  in  overcoming  the  convolutional  ambiguity  in 
the  blind  problem 
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ALGEBRAIC  METHODS  FOR  IMPROVED  BLIND  RESTORATION 
OF  ADAPTIVE  OPTICS  IMAGES  OF  SPACE  OBJECTS 

Brian  D.  Jeffs1 


1  Introduction 

This  document  is  the  final  report  of  progress  made  during  calendar  year  1999  under  the  AFOSR 
Summer  Research  Extension  Program,  Subcontract  No.  98-0813.  Participants  include  Dr.  Brian 
D.  Jeffs  of  Brigham  Young  University  as  principal  investigator,  and  graduate  students  Brent 
Chipman  and  Lisha  Li.  This  work  was  accomplished  In  collaboration  with  Dr.  Robert  Q.  Fugate 
and  Dr.  Julian  C.  Christou  of  the  Starfire  Optical  Range,  Air  Force  Research  Laboratories, 
Kirtland  Air  Force  Base,  New  Mexico. 

The  focus  of  this  research  program  is  algorithm  development  for  blind  image  restoration  of 
multiframe  adaptive  optics  (AO)  telescope  images.  In  particular,  methods  of  primary  interest 
are  those  algorithms  capable  of  removing  residual  blur  from  point-source  images  acquired  using 
the  SOR  ground  based  1.5m  and  3.5m  AO  telescope  systems.  Previous  AFOSR  supported 
work  with  Starfire  Optical  Range  emphasized  blind  restoration  of  AO  images  of  satellites.  The 
restoration  algorithms  used  imaging  models  which  exploited  the  known  man-made  structure  of 
the  true  satellite  image  in  order  to  achieve  higher  resolution  restorations  which  better  recovered 
the  sharp  edges  and  geometric  structures  of  satellites.  The  work  reported  here  extends  these 
ideas  to  the  case  of  blurred  point-like  images,  such  AO  images  of  star  fields  where  separately 
resolving  and  localizing  distinct  stars  which  have  been  blurred  together  is  the  goal.  These  blind 
restorations  are  achieved  by  exploiting  the  known  point-like  image  structure  into  the  restoration 
algorithm,  so  that  point-like  results  are  favored.  Two  new  methods  for  this  class  of  image 
restoration  problems  were  developed,  and  are  reported  here.  The  first  method  yields  a  closed- 

1  Sponsored  by  the  Air  Force  Office  of  Scientific  Research,  Bolling  Air  Force  Base,  DC.  The  author  acknowledges 
the  contributions  to  this  work  from  Dr.  Julian  C.  Christou  and  Dr.  Robert  Fugate,  Starfire  Optical  Range,  Air 
Force  Research  Laboratory.  Also,  Mr.  Brent  Chipman  contributed  much  to  this  work  while  a  graduate  student 
at  Brigham  Young  University. 
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form  algebraic  solution  based  on  subspace  decomposition  techniques.  The  second  method  is 
based  on  Bayesian  estimation  theory  with  a  statistical  model  for  the  sparse  image  field  employed 
as  prior  information. 

2  Subspace  Decomposition  Methods  for  Blind  Point  Source 
Restoration 

In  this  section  we  assume  one  frame  of  blurred  observation  data  is  available  with  no  knowledge 
of  the  blurring  psf.  The  underlying  true  image  is  assumed  to  be  sparse,  consisting  of  a  relatively 
small  number  of  discrete  point-like  objects  (e.g.  stars)  with  unknown  intensity,  and  position.  We 
will  also  assume  the  unknown  blur  is  severe  enough  to  in  some  cases  eliminate  distinct  intensity 
peaks  associated  with  adjacent  points. 

Array  processing  algorithms,  like  ESPRIT  [1]  and  MUSIC,  have  been  used  to  solve  related 
problems  because  they  inherently  deal  with  point  sources.  Recently  the  MUSIC  algorithm  was 
used  to  restore  point-source  images  with  known  blur  [2].  In  this  paper,  we  propose  an  ESPRIT- 
like  algorithm  for  restoration  with  unknown  blur.  The  method  is  non-iterative,  deterministic, 
and  operates  on  a  single  image  frame  with  unknown  blur.  This  formulation  requires  the  fairly 
stringent  constraint  that  the  blurring  psf  be  separable  along  perpendicular  axes.  This  type  of 
blur,  however,  is  not  unrealizable.  Long  exposure  atmospheric  turbulence  blur  has  been  shown 
to  be  approximately  circularly  Gaussian,  which  is  always  separable. 

An  application  of  particular  interest  to  us  is  blind  restoration  of  adaptive  optics  (AO)  telescope 
images  of  star  fields.  The  AO  system  removes  much  of  the  atmospheric-turbulence-induced 
blurring,  but  a  residual  random,  unknown  blur  remains.  The  general  structure  of  the  blur  is 
known  to  be  Lorentzian  [3],  which  is  well  approximated  by  the  sum  of  two  Gaussian  blurs.  With 
this  motivation  we  study  images  blurred  by  non-separable  psf’s  which  can  be  approximated 
with  a  finite  series  expansion  sum  of  separable  blurs.  A  method  is  introduced  for  separately 
estimating  the  x  and  y  positions  of  the  point  sources  in  the  non-separable  case. 
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The  assumed  model  for  the  input  image,  a  point-like  scene  with  d  sources,  is 

d 

f(x,y)  =  -xk)5(y  -yk).  (1) 

k= 1 

The  kth  point  source  has  amplitude  a *  and  position  (rr^y*).  The  observed  image  is  given  by 


g(x,  y)  =  h(x ,  y)  *  f{x ,  y)  +  t?(:e5  y). 


(2) 


Here  the  symbol  V  signifies  two  dimensional  convolution,  h(x,y)  is  the  spatial  domain  psf,  and 
r](x,  y)  is  additive  noise.  The  problem  at  hand  is  to  estimate  point-source  position  pairs,  (x*,  y*), 
and  their  associated  amplitudes. 

2.1  Restoration  with  Separable  Blur 

For  a  separable  blur,  h(x^y)  =  hx(x)hy(y).  Assuming  that  the  output  data  is  available  in  a 
sampled  representation,  y[m,n],  the  two-dimensional  discrete  Fourier  transform  of  the  output 
can  be  written  as  y  (£,£),  where  £  and  (  represent  the  x  and  y  spatial  frequency  variables 
respectively.  Sampling  this  2-D  DFT  with  N  spatial  frequency  samples  along  each  axis,  ^  = 
(i  =  with  i  =  1  •  •  •  AT,  leads  to  a  matrix  representation  of  the  frequency  domain  data: 

G ij  =  g(£jXi)-  Incorporating  equations  (1)  and  2)  into  this  sampled  frequency  domain  model 
yields  the  following  matrix  representation, 


G  =  H2/V(y)AVr(x)H^,  where 
Hy  =  Diag{[hy{(i),--- ,hy{(N)]T}, 
V(y)  =  (v(yi)  •  •  •  v(yd)] , 
v(yfc)  =  •  •  •  e~^NVk^  , 

A  =  Diag{[au---  ,ad]T}. 


(3) 


and  where  Diag{x.)  indicates  a  diagonal  matrix  formed  from  vector  x.  Hx,  V(x),  and  v(a;/t) 
are  defined  similarly.  Note  that  the  kth  element  of  x,  a:*,  is  linked  to  corresponding  element  y 
of  y  in  specifying  x  and  y  positions  of  a  single  source  point. 
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Following  the  development  of  Swindlehurst  and  Gunther  (which  they  used  for  blind  resolution 
of  overlapping  echoes  in  digital  communications)  [4],  we  compute  the  singular  value  decomposi¬ 
tion  G  =  USVF,  where  U  and  V  are  N  xd  and  S  is  d x d.  In  the  noiseless  case,  signal  subspace 
matrices  are  defined  as 

E „  =  US1/2,  and  Ex  =  VS1/2,  (4) 

where  V  =  conj(V).  These  matrices  can  be  related  to  the  model  in  the  following  fashion 

E„  =  HyV(y)A1/2T 

Ex  =  Hx  V  (x)  A1/2  (Tt)  -1 ,  (5) 


where  T  is  some  unknown  full  rank  matrix.  Let  us  define  Vi(x)  as  the  first  N  —  8  rows  of 
V(x)  and  V2(x)  as  the  last  N  -  8  rows  of  V(x).  This  is  equivalent  to  selecting  two  shifted 
sensor  sub-arrays  in  ESPRIT  based  direction  of  arrival  estimation.  Since  the  V  matrices  are 
Vandermonde,  V2(x)  =  Vi(x)3?x  where 

$x  =  Diag  [l,  e~^nSx^N ,  •  ■  • ,  e~j27rSx^Nj  .  (6) 

Similarly  defining  1  and  2  matrices  for  Hx,  H^,  Ex,  Ey,  and  V(y)  leads  to  the  following  relations: 


Ey,  2 

% 

=  T_1$j,T, 

(7) 

EXj2 

=  HX)2EXji^, 

’S'x 

=  T-1$XT. 

(8) 

This  observation  implies  two  equivalent  constraints  on  the  solution.  First,  T,  and  Ty  must 
have  an  equivalent  set  of  eigenvectors  and  second,  \&x.  and  must  commute. 
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2.1.1  The  Optimal  Solution 


An  optimal  solution  to  this  problem  (in  the  least  squares  sense)  can  be  defined  as  the  solution 
to  a  constrained  minimization  problem.  First,  the  following  definitions  are  made 


Ay  —  and  Ax  —  HX)2  *HX)j. 


(9) 


The  solution  is 


Ay,  Ax,  4>y,  ¥x  = 

argminAl(iI)¥y>I 


AyFlyfi 

^^:eEX}2 


Ey,i 

B.,i 


2 


F 


subject  to  the  constraint 


(10) 


eigv{Vy)  =  or  =  0, 


(11) 


where  eig{A}  denotes  the  matrix  of  eigenvectors  of  A.  Once  the  \FX  and  matrices  are  known, 
the  position  estimates  are  given  by 


Xi  = 


Nl\x>i 

2t tS  ’ 


Vi  = 


27:6 


(12) 


where  Ax?*  and  Xy jt-  are  the  ith  eigenvalues  of  and  respectively. 


2.1.2  A  Suboptimal  Solution 

Though  the  solution  presented  above  is  optimal,  it  is  difficult  to  compute  the  joint  minimization 
of  equation  (10).  The  following  suboptimal  approach  performs  the  minimizations  sequentially. 
Ignoring  the  constraint,  it  is  easy  to  solve  for  ^ y  and  Ay  [4], 

8y  =  arg  min  Sjf  [p^  ,  ©  (E„i2E £2)T]  Sy  (13) 
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5X  =  C  lT)<f>x,  and 

(j)x  =  argmin^(E  —  Di?C_1D)</)x  (16) 

<t>x 

where  </>x  =  diag(*x),  Sx  =  diag{ Ax),  E  =  I  0  RHR,  C  =  I  ©  (QQff)r,  D  =  R  ©  Q,  Q 
EX)2(T)t,  and  R  d=  Ex>i(T)t.  It  is  easy  to  solve  for  (f>x  with  a  singular  value  decomposition. 
X{  positions  are  computed  directly  from  (j)x  using  equation  (6),  and  Sx  is  used  to  estimate  H^. 
The  point  source  y  positions  are  similarly  computed  by  switching  the  roles  of  x  and  y  in  the 
derivation  above. 


2.2  Non-Separable  Blurring  Functions 

When  h(x,y)  is  not  separable,  the  above  method  fails,  which  limits  use  of  the  algorithm  to 
a  specific  class  of  images.  This  section  introduces  a  method  that  works  with  non-separable 
blurs.  Suppose  we  have  a  point-source  image  that  is  blurred  by  a  psf  with  exactly  p  separable 
components,  i.e.  h(x,y)  =  hyi{x)hXi{y)-  Note  that  for  sampled  images,  such  a  finite  series 
expansion  of  the  blur  into  separable  components  is  always  possible  by  letting  hyi(m)  and  hXi(n) 
be  the  Fourier  basis  functions.  The  received  image  can  be  modeled  as 

d  v 

9{x->  V )  “  'y  y  ttkhxi{x  —  %k)h,yi{y  Vh)  (17) 

fc=l  i= 1 
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Taking  a  1-D  DFT  of  the  columns  and  taking  N  samples  as  before  yields 


G,  =  EH*V(y)A*  +  NF.  (18) 

1  =  1 

Each  of  the  separable  components  are  of  rank  d ,  so  the  total  rank  in  the  noise  free  Gy  must  be 
q  <  pd.  Letting  E  be  the  first  q  singular  vectors  from  the  SVD  of  Gy  leads  to 

[H„iV(y)|  •  •  •  |HyPV(y)]  €  span{E}.  (19) 


Thus  for  each  i  there  exists  a  q  x  d  matrix  Tj  such  that 


ETj  =  HyjV(y). 


(20) 


This  can  be  recast  into  the  least  squares  problem 


H yi,  T i,  y  =  arg  min  ||Hj,iV(y)  -  ETj||^ . 


(21) 


Using  the  sub-optimal  sequential  estimation  approach  to  solve  for  T *  first  gives 


Ti  =  (EEW)  1  EtfHytV(y). 


(22) 


Reinserting  this  solution  into  equation  (21)  yields 


Hyi,y  =  arg  min 
H,y 


P^HyiV(y)||F, 


(23) 


which  can  be  expanded  in  exactly  the  same  manner  as  the  previous  section’s  solution  for  6  to 
obtain 


hyi,y  =  argmaxh^0(y)h!/i, 


h,y 


where 


®(y)  =  (eEh)  ©  (v(y)Vff  (y))T . 


(24) 


(25) 
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If  the  vectors  h yi  are  constrained  to  be  of  unit  magnitude,  the  maximum  value  of  the  quadratic 
form  in  equation  (24)  will  be  Amax,  the  maximum  eigenvalue  of  ©(y).  Thus,  regardless  of  the 
values  of  the  true  h yi  vectors,  the  solution  for  positions  can  be  expressed  independent  of  the 
unknown  blur  with  the  following  reformulation 

y  =  arg  max  {max  eigval  (®(y))}  .  (26) 

y 

This  appears  to  be  a  daunting  equation  to  solve  numerically.  At  first  glance  it  seems  as  though 
there  is  no  hope  without  an  exhaustive  search  over  the  d  dimensional  parameter  space  in  y,  (i.e. 
compute  0(y)  for  all  possible  configurations  of  d  points  positioned  along  the  y  axis,  calculating 
eigenvalues  at  each  step.)  It  is  possible  however  to  solve  the  problem  with  at  most  2  scans  over 
a  one  dimensional  parameter  space  by  exploiting  underlying  structure. 

While  equation  (20)  holds  for  the  vector  y  consisting  of  all  d  of  the  y  positions  for  point 
sources,  it  also  holds  for  a  shorter  y  containing  at  least  two  elements  (points).  Correspondingly, 
equation  (26)  will  achieve  maxima  for  any  size  vector  y  as  long  as  all  the  specified  positions 
match  some  subset  of  points  in  truth  image  f(x,  y).  In  particular,  if  only  two  points  are  chosen, 
the  largest  eigenvalue  of  0(y)  will  achieve  a  local  maximum  whenever  the  two  y  parameters 
are  separated  by  any  distance  which  exists  between  some  actual  point-source  pair.  We  note 
that  in  any  truly  blind  restoration,  only  relative  distance  between  features  can  be  recovered. 
Absolute  position  is  irrelevant  because  global  shifts  with  respect  to  a  reference  frame  can  also 
be  interpreted  as  arising  from  a  non-centralized  blur  psf. 

For  example,  consider  an  image  with  5  points  located  at  positions  y  =  {10,20,25,37,39}. 
The  set  of  all  pairwise  differences  in  position  is  {2, 5, 10, 12, 14, 15, 17, 19,  27, 29}.  Defining  a  two 
point  position  vector  as  y  =  [0 ,  y]r,  turns  equation  (26)  into  the  following  1-D  optimization 
problem  in  y 

f{y )  =  max  eigval(Q(y ) ) .  (27) 

The  plot  of  /  as  a  function  of  y  is  called  the  difference  scan ,  which  will  manifest  peaks  at  each 
legitimate  difference.  It  is  possible  to  reduce  the  computation  by  appealing  to  a  lemma  found  in 
[4]  whereby  we  can  reduce  the  matrix  0(y)  to  a  size  of  q  x  g,  which  is  smaller  than  an  N  x  N, 
Armed  with  the  difference  scan,  we  can  obtain  a  solution  for  the  position  estimates.  The  two 
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points  that  are  separated  by  the  maximum  difference  are  uniquely  determined.  To  prove  this 
property,  assume  there  are  two  sets  of  point  pairs,  (yi,y2)  and  (y3,y4),  which  are  separated  by 
the  same  (maximum)  distance,  Dmax.  For  this  condition 

V2  ~  yi  =  D  max’)  and  y4  -  y3  =  D 

max  *  (28) 

If  2/1  <  2/3  <  V2  then  clearly  y4  —  yi  >  Dmax  which  violates  the  premise  that  Dmax  was  the 
maximum  difference.  The  other  possible  relative  position  cases  lead  to  a  similar  conclusion. 
Therefore,  the  local  peak  in  equation  (27)  corresponding  to  the  largest  y  value  localizes  the 
unique  point  pair  distance  Dmax.  With  this  first  pair  separation  uniquely  determined,  each 
other  position  can  be  uniquely  determined  between  them.  Now,  define  the  single  parameter 
point  triple 

y  =  [Oilman  2/]  >  (29) 

and  the  solution  spectrum,  /(y),  will  have  unique  peaks  at  all  correct  positions,  including  0  and 
Dmax,  when  scanning  the  single  scalar  parameter  y. 

2.3  Results  for  the  Subspace  Approach 

The  first  example  shown  in  Figure  1  is  for  the  separable  blur  method.  The  observed  image 
contains  six  stars  convolved  with  a  circularly  symmetric  Gaussian  blurring  function  with  a 
standard  deviation  of  3.5  pixels.  Two  sets  of  stars  are  very  much  blurred  together.  The  top 
left  frame  shows  the  blurred  data  with  white  Gaussian  noise  added  at  a  level  of  40  dB  peak 
SNR.  The  top  right  frame  shows  the  position  estimates  as  asterisks,  and  the  correct  positions 
as  diamonds.  These  results  are  obtained  with  the  y  positions  estimated  without  regard  to  the 
constraint  of  equation  (11).  The  bottom  left  image  shows  the  converse  case,  where  estimates  of 
the  x  positions  do  not  include  the  constraint.  Note  that  position  estimates  are  quite  accurate 
(average  error  less  than  two  pixels),  which  is  remarkable  given  the  noise  level  and  complete  lack 
of  knowledge  about  the  blur,  h(x,y),  other  than  that  it  is  separable. 

The  second  example  shown  in  Figure  2  illustrates  the  non-separable  method.  The  same  input 
image  is  used,  except  this  time  the  blur  is  a  non-separable  Lorentzian  function.  The  top  left 
frame  shows  the  blurred  image,  while  the  top  right  shows  the  input  image  itself.  Note  that  the 
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this  truth  image  corresponds  to  a  smaller  window  centered  on  the  blurred  points  in  the  observed 
image,  and  is  thus  presented  at  a  lower  resolution  scale.  The  bottom  two  graphs  are  the  solution 
spectra  for  the  y  and  x  positions  respectively.  We  see  that  the  peaks  of  the  spectra  indeed 
correspond  to  the  actual  positions  represented  in  the  input  image. 


Figure  1:  Six  star  example  at  PSNR=40  dB.  Top  left:  Blurred,  noisy  output.  Top  right: 
Position  estimates  with  y  fixed.  Bottom  left:  Position  estimates  with  x  fixed.  Bottom  right: 
Blur  estimate. 


3  Bayesian  Methods  for  Blind  Point  Source 
Image  Restoration 

In  this  section  we  present  a  Bayesian  estimation  theory  based  approach  for  solving  the  blind 
multiframe  point-source  image  restoration  problem.  We  assume  one  or  more  frames  of  blurred 
observation  data  are  available  with  no  knowledge  of  the  blurring  psf,  and  that  blur  psf’s  may 
be  different  from  frame  to  frame  if  multiple  observations  are  available 
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Figure  2:  Six  star  example  200  dB.  Top  left:  Blurred,  noisy  image.  Top  right:  Truth  image. 
Bottom  left:  y  positions  solution  spectrum.  Bottom  right:  x  position  solution  spectrum. 

An  application  of  particular  interest  to  us  is  blind  restoration  of  adaptive  optics  (AO)  telescope 
image  sequences  of  star  fields.  The  AO  system  removes  much  of  the  atmospheric  turbulence 
induced  blurring,  but  a  residual  random,  unknown  blur  remains  that  changes  from  frame  to 
frame  in  an  image  sequence  over  a  period  of  milliseconds.  Though  the  gross  structure  of  the  blur 
is  known  on  average  [6],  the  specific  form  of  each  individual  blur  cannot  be  easily  ascertained,  and 
is  thus  best  modeled  with  the  MRF  approach,  proposed  in  this  paper.  Identifying  individual  stars 
in  dense  star  clusters,  forming  accurate  photometry  estimates,  and  computing  star-positions  to 
sub-pixel  accuracy  are  primary  goals.  For  example,  precise  measurement  of  relative  positions 
can  help  identify  the  “wobble”  associated  with  stars  orbited  by  massive  planets.  This  situation 
lends  itself  well  to  the  blind  restoration  technique  presented  in  this  paper. 

Bayesian  maximum  a  posterior  (MAP)  estimation  has  been  shown  to  be  effective  in  blind 
restoration.  In  particular,  Jeffs,  Hong,  and  Christou  [7]  have  recently  demonstrated  the  effec¬ 
tiveness  of  generalized  Gauss  Markov  random  fields  (GGMRF)  in  blind  restoration  of  extended 
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objects.  Here  both  the  source  and  blur  were  modeled  as  GGMRF’s  which  have  a  parametric 
form  allowing  a  great  variety  of  image  representations,  including  hard  edged  fields  typical  of 
real  images  and  smooth  fields  typical  of  blurring  point  spread  functions.  However,  the  GGMRF 
model  is  not  well  suited  to  point-like  sparse  images.  A  Markov  random  field  model  which  favors 
sparse  solutions  is  essential  if  high  resolution  restorations  and  accurate  point  localizations  are  to 
be  achieved,  particularly  in  the  blind  case.  The  ability  to  exploit  known  structure  in  the  problem 
and  impose  a  sparse  form  on  the  solution  is  essential  in  overcoming  convolutional  ambiguity  in 
the  blind  problem.  Blind  restoration  is  a  highly  ill-posed  inverse  problem,  and  algorithms  which 
incorporate  known  image  structure  in  solutions  will  invariably  perform  better. 

Phillips  and  Leahy  [8],  have  presented  an  MRF  model  that  exploits  the  sparse  nature  of  point 
source  input  images  in  the  context  of  MEG-based  imaging.  Their  model  involves  a  dual  field 
representation:  first,  a  binary  activity  process  determines  which  pixels  have  non-zero  amplitudes, 
then  a  Gaussian  amplitude  process  represents  active  point  intensity  levels.  We  demonstrate  that 
this  model  can  be  effectively  extended  to  the  blind  point  source  restoration  case.  In  fact,  the 
prior  information  provided  by  this  sparse  model  image  prior  pdf  is  the  key  to  overcoming  inherent 
ambiguity  when  blur  psfs  are  unknown. 

3.1  Problem  Formulation 

We  adopt  the  following  image  observation  model  for  both  single  and  multiple  frame  data  repre¬ 
sentation 


I  =  m  +  fj  (30) 

U  = 

where  M  is  the  number  of  frames,  gi,  f,  and  rji  are  vectors  formed  by  column  scanning  the 
2-D  images  of  the  ith  observation  frame,  the  true  image,  and  the  ith  noise  frame  respectively. 

is  the  doubly  block  Toeplitz  convolution  matrix  formed  from  the  ith  frame  psf,  hi.  g  is 
the  extended  vector  formed  by  concatenating  all  M  column  scanned  observation  frames,  and 
K  represents  all  M  distinct  frame  blur  psf’s  as  a  single  system  matrix.  This  formulation  also 
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works  for  the  single  frame  case  by  setting  M  =  1. 

Assuming  f  and  H  are  statistically  independent,  the  blind  MAP  restoration  problem  may  be 
stated  as 

=  argmaxpg|M(g|f,^)p/(f)p/l('H)  (31) 

I 

We  will  assume  that  ^(77)  is  zero  mean,  i.i.d.  Gaussian.  This  implies  that  pg\fth( g|f,^)  is 
i.i.d.  Gaussian  with  a  mean  of  Hf . 

In  order  to  model  f  as  a  sparse  MRF,  we  follow  the  development  of  Phillips  and  Leahy  by 
defining 


f  =  Xz  (32) 

where  X  is  a  diagonal  matrix  with  elements  of  either  0  or  1.  The  vector  z  is  vector  of  amplitudes. 
The  vector  x  =  diag{X },  represents  an  indicator  process  that  determines  whether  or  not  a 
particular  pixel  is  active.  In  the  solution  of  equation  (31)  we  make  the  replacement 

Pf{  f)  =  z)  =  Px{x)Pz{  z)  (33) 

where  we  assume  the  indicator  process  and  the  amplitude  process  are  independent.  The  indicator 
function  can  be  modeled  as  a  binary  Markov  random  field  whose  probability  density  function 
follows  a  Gibb’s  distribution, 


PW  =  ^exp(-V(x)))i,G{0,l}  (34) 

where  if  is  a  normalizing  constant  and  the  Gibbs  distribution  potential  function,  F(x),  is  given 

by 


V(x)  =  Y^diXi  +  fkCi{xi,Xj,j  eAfi}  (35) 

i 

where  and  $  are  weighting  constants  and  C{  is  a  clustering  function  that  operates  on  pixels 
in  the  neighborhood  Mi  [8].  Arguing  that  there  is  no  reason  to  suspect  any  clustering  a  priori 
of  adjacent  pixels  in  star  images  (clusters  may  exist,  but  individual  stars  would  be  separated  by 
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Figure  3:  Image  to  be  restored.  Top  left:  Actual  blur  psf  formed  as  an  elliptical  Lorentzian 
function  on  a  GGMRF  residual  halo  producing  a  realistic  AO  low-pass  “mottled”  halo  field. 
Top  right:  Circularly* symmetric  Lorentzian  blur  model  to  be  used  as  the  reference  mean,  Hh- 
Bottom  left:  Actual  truth  image.  Bottom  right:  Observed  blurred,  noisy  output  data. 

some  black  space,  which  is  different  from  the  Phillips-Leahy  clustering  model),  we  have  omitted 
the  clustering  term  and  focused  on  the  first  term  in  equation  (35)  which  enforces  sparseness  in 
the  final  image. 

The  density  function  for  the  amplitudes  is  assumed  to  be  Gaussian  in  the  Phillips-Leahy 
approach.  They  are  dealing  with  MEG  processing  in  which  the  amplitudes  can  be  either  positive 
or  negative  and  thus  introduce  a  zero  mean  Gaussian.  In  star  images,  where  we  deal  only  with 
intensities,  negative  amplitudes  are  unacceptable.  This  can  be  enforced  by  adding  a  mean  to  the 
Gaussian,  mz.  Another  method  is  to  simply  draw  the  amplitudes  from  a  uniform  distribution 
over  some  positive  range  and  let  the  noise  term  in  Equation  (31)  dictate  the  final  values. 


Figure  4:  Restoration  results.  Top  left:  Estimated  blur.  Top  right:  Restored  activity  matrix. 
Bottom  left:  Restored  input  image.  Bottom  right:  Restored  image  blurred  with  restored  blur 
psf.  Compare  with  Figure  3. 


The  blur  pdf  is  modeled  as  a  GGMRF  with  density 


Ph{%)  ^XP  |  ^  ^  1  Ph,s  | ^ 


01  ^  cs,t|(^s  M h,s )  {ht 

<s,t>ech 


(36) 


where  Ch  is  the  set  of  all  cliques  for  the  blur  neighborhood  system,  q  is  the  blur  GGMRF  shape 
parameter,  5/,  is  the'set  of  all  points  in  the  blur  lattice  over  all  frames  and  cSit  and  ds  are 
neighborhood  influence  weights. 

It  has  been  shown  in  previous  publications,  for  example  [7],  that  q  >  2  does  a  good  job 
modeling  smooth  features  such  as  those  found  in  a  typical  blur  psf.  The  mean,  /xh>g,  allows  the 
restoration  to  maintain  fidelity  with  prior  knowledge  about  the  blur.  For  example,  in  astronom- 


ical  imaging  isolated  stars  in  a  nearby  field  can  be  averaged  to  give  a  reference  mean. 

Combining  the  density  functions  we  can  re-express  equation  (31)  by  taking  the  logarithm  of 
both  sides  and  dropping  the  additive  constants 

M 

X,z,?i  =  arg  min  y  ||gj  -  HjXz||2+  (37) 

x'z'ni= i 

7]C«j  +  A||z-jiz||2 

3 

~  \q  + 

sesh 

a  cs,t\(hs  -  Hs,t)  ~  {ht  -  fJLh,t) \q 

<s,t>ech 

Here  A, 7, a  control  the  relative  influence  as  regularizing  terms  that  the  activity  matrix,  ampli¬ 
tudes,  and  blur  have  on  the  solution. 

Equation  (37)  represents  a  very  complicated  nonlinear  minimization  problem.  The  approach 
taken  to  solve  it  is  simulated  annealing  [9,  10],  specifically  the  Metropolis  algorithm  [11].  A, 7, a 
are  set  manually  and  adjusted  for  best  restoration  performance. 

3.2  Results  for  Bayesian  Blind  Point  Source  Restoration 

In  this  section  we  present  examples  of  the  new  blind  algorithm  using  simulated  adaptive  optics 
telescope  data.  In  Figure  (3),  the  actual  data  for  the  first  example  is  presented.  This  shows  a 
blur  with  region  of  support  that  is  15  x  15  pixels.  The  blur  is  a  rotated  Lorentzian  function 
with  an  elongated  axis.  This  has  been  shown  to  be  an  excellent  model  for  AO  residual  blur  [6]. 
The  reference  mean  used  for  the  blur  is  a  circularly  symmetric  Lorentzian  shape,  with  different 
radius  than  either  axis  of  the  actual  blur.  The  truth  image  is  also  shown  consisting  of  ten 
isolated  points.  The  bottom  right  image  shows  the  resulting  blurred,  noise  corrupted  output. 
The  noise  is  zero-mean  white  Gaussian  noise  at  a  level  of  32  dB  peak  SNR. 

Note  that  only  one  observed  frame  is  used  in  this  restoration  example.  This  is  actually  a  more 
difficult  problem  than  the  case  of  multiple  frames  with  different  blurs  for  each  frame.  In  the 
multiframe  case,  the  diversity  provided  by  distinct,  unknown  blurs  aids  in  the  blind  restoration 
because  only  the  true  image  is  common  to  all  frames. 

Figure  (4)  shows  the  results  of  the  algorithm  described  above.  The  blur  has  been  estimated 
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quite  accurately.  The  main  axis  is  quite  clearly  defined  in  the  restoration  and  the  extent  of 
the  restored  blur  matches  the  actual  blur  very  well.  The  source  restoration  matches  the  actual 
image  to  a  high  degree.  However,  two  points  in  the  L-shaped  structure  of  the  original  source 
were  blurred  into  a  single  point  in  the  solution  and  some  positions  are  inaccurate  by  a  single 
pixel. 


Figure  5:  Second  example.  Top  left:  Observed  blurred,  noisy  output  data.  Top  right:  Actual 
truth  image.  Bottom  left:  Actual  blur  psf.  Bottom  right:  Circularly  symmetric  Lorentzian 
reference  mean. 


The  actual  data  for  the  second  example  is  found  in  Figure  (5).  Here  we  have  a  similar  blur 
with  a  input  image.  Figures  (6)  and  (7)  show  two  different  restorations  on  the  data.  In  each 
of  the  restorations,  the  blur  is  only  slightly  elongated  suggesting  too  much  weight  on  the  blur 
self  term  (i.e.  on  /i^).  Each  of  the  input  image  restorations  are  fairly  good.  The  first  misses  a 
point  to  the  bottom  right  of  the  input  image  and  some  other  positions  are  off  by  about  a  pixel. 
The  second  get  all  points  but  in  the  upper  portion  of  the  input  image,  one  of  the  points  is  split. 
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Figure  6:  Second  example,  first  restoration.  Top  left:  Restored  image  blurred  with  restored  blur 
psf.  Top  right:  Restored  activity  matrix.  Bottom  left:  Restored  input  image.  Bottom  right: 
Restored  blur  psf. 

This  is  likely  due  to  very  circular  restoration  of  the  blur.  These  problems  would  likely  be  solved 
with  better  weights  oxy  the  various  priors. 

4  Conclusions 

This  report  has  detailed  the  progress  made  under  support  from  1999  AFOSR  Summer  Research 
Extension  Program  for  research  in  multiframe  blind  restoration  of  space  objects  images  from 
adaptive  optics  systems.  This  effort  has  been  very  fruitful,  and  some  very  promising  new  tech¬ 
nologies  have  been  identified.  It  is  expected  that  this  research  effort  will  lead  to  higher  resolution 
restorations  of  AO  images  with  less  computational  burden  than  other  blind  algorithms. 

We  have  presented  two  methods  of  blind  point-source  image  restoration  based  on  determin- 


Figure  7:  Second  example,  second  restoration.  Top  left:  Restored  image  blurred  with  restored 
blur  psf.  Top  right:  ^Restored  activity  matrix.  Bottom  left:  Restored  input  image.  Bottom 
right:  Restored  blur  psf. 

istic,  subspace  decomposition  algorithms,  and  a  third  method  using  sparse  statistical  models 
as  the  basis  of  a  Bayesian  estimation  approach.  For  the  subspace  approach,  the  underlying 
point-source  image  is  recovered  from  a  single  observation  frame  without  any  prior  knowledge 
about  the  blurring  function  (in  the  case  of  the  non-separable  algorithm).  If  it  is  known  that 
the  blur  psf  is  separable,  then  an  even  more  computationally  efficient  algorithm  is  available. 
Other  known  point  source  restoration  methods  are  iterative,  or  require  knowledge  of  the  blur,  or 
require  multiple  observation  frames  with  distinct  blur,'  or  are  based  on  statistical  signal  models, 
or  involve  lengthy  simulated  annealing  optimization  codes  [2,  5,  12,  13].  We  believe  the  new 
methods  presented  here  are  superior,  first  in  elegance  and  second  in  speed.  We  know  of  no 
other  deterministic  blind  method  which  is  optimized  for  point-source  images.  The  new  methods 


involve  extension  and  adaptation  of  digital  communications  algorithms,  originally  developed  for 
resolving  overlapping  multipath  echoes,  so  that  they  can  be  used  in  2-D  blind  imaging  prob¬ 
lems.  Future  work  should  include  further  attention  to  the  issue  of  linking  the  separate  x  and 
y  spectrum  scanned  position  estimates  in  the  non-separable  algorithm.  Methods  of  extracting 
blur  psf  and  point  amplitude  estimates  from  the  non-separable  algorithm  are  also  under  study. 

The  Bayesian  approach  was  also  shown  to  be  a  viable  tool  in  solving  the  blind  point  source 
image  restoration  problem.  The  restoration  examples  shown  above  demonstrate  the  power  of 
the  model  described  in  this  report  to  emphasize  the  sparse  character  of  a  source  image.  Though 
GGMRFs  have  been  shown  to  work  well  for  restoring  extended  objects,  the  model  does  not 
allow  the  user  to  explicitly  enforce  point-like  structure  on  the  solution.  The  specific  point- 
source  prior  of  Phillips  and  Leahy  in  conjunction  with  the  GGMRF  prior  for  the  psf  leads  to 
blind  restoration  solutions  which  are  truly  sparse  and  are  an  excellent  estimate  of  the  truth.  We 
have  demonstrated  that  the  Phillips  -  Leahy  MRF  model  for  point  sources  is  well  suited,  with 
minor  modifications,  to  blind  star  field  image  restoration.  Previously  this  model  had  been  used 
only  with  known  blurring  or  system  functions. 

4.1  Additional  Accomplishments  Supported  by  the  SREP  Grant 

This  section  summarizes  accomplishments  resulting  from  the  AFOSR  Summer  Faculty  Research 
Extension  Program  support  of  my  research  collaboration  with  Starfire  Optical  Range  at  Kirtland 
Air  Force  Base.  These  activities  are  in  addition  to  the  research  progress  reported  above,  and 
include  work  completed  from  the  end  of  1998  until  the  present. 

List  of  Additional  Accomplishments: 

•  Technical  Papers  Resulting  From  SOR  Work: 

1.  Brent  A.  Chipman  and  Brian  D.  Jeffs,  “Blind  Multiframe  Point  Source  Image  Restora¬ 
tion  Using  MAP  Estimation, 51  Conference  Record,  Thirty-Third  Asilomar  Conference 
on  Signals,  Systems,  and  Computers ,  Nov.  1999,  Pacific  Grove  CA. 

2.  Brent  A.  Chipman  and  Brian  D.  Jeffs,  “Blind  Point-Source  Image  Restoration  Using 
Subspace  Techniques,”  to  appear  in  Proceedings  of  IEEE  International  Conference 
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on  Acoustics,  Speech,  and  Signal  Processing,  ICASSP-2000,  Istanbul  Turkey,  June 

2000. 

•  Master’s  thesis  resulting  from  this  funded  research:  Brent  Chipman,  Methods  of  Blind 
Point  Source  Restoration,  Brigham  Young  University,  Dec.  1999. 

•  Student  support:  This  research  program  has  had  a  profound  effect  on  two  of  my  graduate 
students.  Brent  Chipman  has  recently  graduated  with  a  Master’s  degree,  and  a  fine  thesis 
reporting  his  work  in  AO  image  post  processing.  Also,  Lisha  Li  has  recently  started  Ph.D. 
dissertation  work  based  on  this  research  program. 

•  Ongoing  work:  includes  a  collaborative  research  effort  with  Julian  Christou  and  Stuart 
Jefferies  of  the  University  of  Arizona.  We  are  currently  preparing  a  proposal  for  AFOSR 
for  continued  work  in  advanced  algorithms  for  blind  restoration  of  AO  images  of  satellites. 

I  expect  that  my  research  in  this  fruitful  area  and  collaboration  with  Starfire  Optical  Range 
will  continue  for  a  number  of  years. 
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Abstract 

A  model  of  a  piezoceramic  transducer  coupled  to  a  cylindrical  cavity  is  derived  for  the  purpose 
of  studying  active  damping  of  interior  noise.  The  model  is  nondimensionalized  to  determine 
the  critical  parameters  that  affect  the  structural-acoustic  coupling  of  the  transducer.  The 
nondimensional  model  yields  four  critical  parameters:  the  mass  ratio  of  the  acoustic  piston 
to  the  beam  flexure;  the  actuator-to-flexure  stiffness  ratio;  the  acoustic-to-flexure  stiffness 
ratio;  and  the  frequency  response  of  the  acoustic  load.  Nondimensional  expressions  for  these 
parameters  are  derived  as  a  means  of  relating  the  transducer  design  to  the  structural-acoustic 
response.  Numerical  analysis  indicates  that  minimizing  the  mass  ratio  and  making  the  stiff¬ 
ness  ratios  approximately  unity  results  in  sufficient  structural-acoustic  coupling.  Simulations 
of  velocity  feedback  demonstrate  the  feasibility  of  the  concept.  Applications  for  this  work 
include  the  development  of  lightweight  transducers  for  problems  such  as  interior  noise  in 
launch  vehicle  payload  fairings. 
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Introduction 


Interior  noise  control  is  important  for  a  number  of  aerospace  applications.  Previous  research 
has  focused  on  understanding  the  fundamental  mechanisms  of  sound  transmission  into  en¬ 
closed  spaces  [1,  2].  Experimental  results  have  been  presented  for  simple  structures  such  as 
rectangular  enclosures  and  cylinders,  and  for  more  complicated  geometries  that  model  the 
structural-acoustics  of  aircraft. 

The  motivation  for  this  work  is  interior  noise  control  for  launch  vehicle  payload  fairings. 
One  of  the  critical  features  of  this  aerospace  application  is  the  need  for  lightweight  trans¬ 
ducers  for  active  control  of  sound.  Although  weight  is  a  primary  design  parameter  for  all 
applications,  it  is  especially  important  in  launch  vehicle  noise  control  because  of  the  cost 
associated  with  placing  payloads  into  orbit.  Current  cost  estimates  for  launch  vehicles  are 
approximately  $10,000  per  pound  of  payload  capacity.  Weight  added  to  the  vehicle  for  the 
purpose  of  active  noise  control  will  reduce  the  payload  capacity  of  the  vehicle.  Therefore  the 
feasibility  of  active  noise  control  for  payload  fairings  is  strongly  coupled  to  the  development 
of  lightweight  transducers  for  controlling  the  interior  noise. 

Most  of  the  experimental  work  in  interior  noise  control  has  utilized  speakers  as  the  acous¬ 
tic  control  source.  Speakers  are  a  mature  technology  that  provide  the  necessary  stroke  and 
acoustic  power  for  noise  control  applications.  A  conventional  speaker  utilizes  an  electromag¬ 
netic  actuator  as  the  driving  element,  a  cone  as  the  acoustic  element,  and  some  measure 
of  stiffness  for  mechanical  stability.  The  moving  mass  of  the  speaker  and  the  mechanical 
stiffness  create  a  mechanical  resonance  that  varies  depending  on  the  speaker  design.  Typical 
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values  for  this  resonance  are  20-50  Hz  for  bass  speakers  and  200-300  Hz  for  smaller  designs. 

Although  speakers  offer  many  advantages  for  noise  control,  they  also  have  drawbacks  that 
are  inherent  in  the  speaker  design.  The  mechanical  resonance  of  the  speaker  introduces  actu¬ 
ator  dynamics  that  can  be  problematic  for  feedback  control.  The  phase  lag  associated  with 
the  resonance  can  cause  stability  problems  and  reduce  achievable  performance.  Eliminating 
problems  due  to  mechanical  resonance  was  the  basis  for  the  work  of  Clark  and  Lane  in  the 
development  of  a  ‘constant  volume  velocity  source’  [3]  They  developed  a  technique  for  esti¬ 
mating  the  velocity  of  an  electromagnetic  speaker  for  the  purpose  of  reducing  the  phase  lag 
associated  with  mechanical  resonance.  They  applied  their  technique  to  the  control  of  ducts 
and  the  fuselage  of  a  small  aircraft.  Another  drawback  of  conventional  speaker  technology 
is  the  sheer  weight  of  the  transducer.  The  magnets  required  for  the  the  driving  element  are 
the  heaviest  component  of  the  device.  This  is  particularly  problematic  for  low-frequency 
applications  that  require  large  strokes. 

The  purpose  of  this  work  is  to  investigate  the  use  of  piezoceramic  actuators  as  driving 
elements  for  noise  control  transducers.  The  objective  is  to  develop  a  lightweight  transducer 
that  overcomes  problems  due  to  mechanical  resonance.  Sound  absorption  will  be  accom¬ 
plished  with  feedback  control  utilizing  the  piezoceramic  transducers.  The  control  algorithms 
will  be  simple  and  robust  to  maximize  reliability. 

The  rationale  for  this  work  is  a  recent  paper  by  Leo  and  Limpert  on  the  development  of 
a  self-sensing  technique  for  acoustic  attenuation  [4].  They  developed  and  implemented  an 
algorithm  that  estimated  the  back  pressure  of  an  electromagnetic  speaker  using  a  measure 
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of  the  electrical  impedance.  The  technique  was  effective  in  frequency  ranges  in  which  the 
electrical  impedance  was  coupled  to  the  acoustic  resonances  of  the  enclosed  space.  Outside 
this  range  the  algorithm  was  overly  sensitive  to  small  errors  in  the  implementation  of  the 
estimated  electromechanical  impedance. 

This  work  highlighted  the  importance  of  coupling  the  transducer  to  the  acoustic  field. 
Subsequent  analyses  indicated  that  if  the  transducer  was  coupled  to  the  acoustic  field  then 
there  was  no  compelling  reason  to  perform  self-sensing;  it  was  sufficient  to  utilize  the  trans¬ 
ducer  as  a  collocated  control  element  for  active  damping. 

Therefore  this  work  will  concentrate  on  identifying  the  critical  parameters  for  developing 
an  acoustic  source  utilizing  piezoceramic  transducers  for  acoustic  feedback  control.  The  next 
section  will  discuss  the  fundamental  concept  and  a  simplified  model  will  be  developed  in  the 
subsequent  section.  A  discussion  of  the  critical  design  parameters  will  be  presented  and 
several  numerical  analyses  will  be  performed  to  demonstrate  the  coupling  of  the  transducer 
to  the  acoustic  field. 


Transducer  Modeling 

Figure  la  is  a  schematic  of  the  transducer  concept.  The  transducer  consists  of  an  acoustic 
piston  and  a  housing  which  is  similar  to  a  conventional  speaker.  The  piston  is  angled  to 
maximize  the  stiffness  per  unit  area  of  the  acoustic  source.  The  actuation  elements  are  a 
set  of  piezoceramic  bimorphs  attached  to  the  bottom  of  the  piston  and  cantileverd  on  the 
housing.  A  single  actuator  and  bimorph  is  called  a  flexure.  The  set  of  flexures  provide  the 
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(a)  (b) 


Figure  1:  Transducer  concept  illustrating  the  flexure  and  actuation  mechanisms:  (a)  top 
view,  (b)  side  view  illustrating  the  assumptions  of  pure  translation. 

mechanical  stiffness  of  the  acoustic  source. 

The  goal  of  this  work  is  to  investigate  the  critical  parameters  of  the  transducer  design. 
To  this  end,  we  will  model  the  structural-acoustic  coupling  of  a  single  flexure  for  a  particular 
acoustic  load.  This  will  allow  us  to  develop  a  set  of  nondimensional  parameters  that  couple 
the  electromechanical  parameters  of  the  flexures  to  the  resonances  of  the  acoustic  cavity. 

Figure  lb  is  a  simplified  model  of  the  transducer  flexure.  For  this  work  we  assume  that 
each  flexure  is  cantilevered  to  the  housing.  Thus  the  displacement  and  rotation  is  zero  where 
the  flexure  is  attached  to  the  housing.  Another  assumption  is  that  the  cone  undergoes  pure 
translation  and  no  rotation.  Thus  the  rotation  of  each  flexure  is  zero  where  the  flexure  is 
attached  to  the  cone.  The  mass  of  the  cone  is  modeled  as  a  concentrated  mass  at  the  end  of 
each  flexure. 
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The  simplified  transducer  model  will  be  separated  into  three  parts:  a  beam  model  for 
the  flexure,  a  one-dimensional  strain  model  for  the  sensors  and  actuators,  and  an  impedance 
model  for  the  acoustic  load.  For  this  work,  the  acoustic  loads  is  assumed  to  be  a  cylindrical 
tube  with  a  rigid  termination.  The  beam  model  will  be  developed  through  an  energy  formu¬ 
lation  and  the  actuator  and  acoustic  loads  will  be  modeled  as  nonconservative  work  terms 
acting  on  the  flexure. 

Flexure  Model 

Each  flexure  is  modeled  independently  using  the  assumed  modes  method.  The  flexure  is 
assumed  to  be  a  thin,  flexible  beam  that  satisfies  the  Euler-Bernoulli  assumptions.  The 
transverse  displacement  of  the  beam  is  denoted  w(x,t).  The  kinetic  and  potential  energy 
terms  for  the  beam  are 


where  p,  A,  E,  and  I  are  the  density,  cross  sectional  area,  modulus,  and  cross  sectional  inertia 
of  the  beam,  respectively.  The  concentrated  mass  mp  models  the  mass  of  the  rigid  piston. 
Piston  inertia  and  the  static  load  due  to  the  piston  are  neglected.  The  nonconservative  work 
due  to  piezoceramic  actuation  is  expressed  as 

w„c  =  m„(,)  - rftJyWM)  (2) 

where  ma  is  the  moment  applied  by  the  actuator,  p  is  the  uniform  back  pressure  applied 
to  the  piston,  and  Ap  is  the  piston  area.  Assuming  that  the  displacement  w  is  separable  in 
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space  and  time  and  expressed  as  a  linear  combination  of  admissible  functions  <f>(x): 


N 

w(x,t)  =  Y/<t>i{x)qi(t)  (3) 

i=l 

where  the  admissible  functions  satisfy  the  geometric  boundary  conditions  of  the  structure, 
we  can  write  the  kinetic  energy  expression  in  matrix  form 

T  =  eYiq(t)T  Mbq(t)  +  \q{t)T  Mcq{t) 

2  (4) 

V  =  siu±-q(tyKiq(t) 

where, 

Mbll  =  MiMt  K 


mp„  =  M«c)0>(fe)  (5) 

K.  fl 

**■ bij  JO  d£, 2  dp 

and  £  =  x/L  is  a  nondimensional  length  and  p  =  The  nonconservative  work  terms  can 
also  be  written  as  a  matrix  expression 


Wnc  _  r^)_BTq(^j  —  p(t)ApBpq(t.) 


(6) 


where 


and, 


Ba 


2)  _  d4>i(f  1)  dfM&l  _  d^Af(6) 


BP  —  [0i  (!)•••  0at(  1)]T 


The  equations  of  motion  are  derived  using  the  Lagrangian  formulation: 

El  1 

pAL  {Mb  +  pMp)  q  +  -jjKbq  =  -Bama  -  ApBpp. 


(7) 

(8) 

(9) 
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Figure  2:  Geometry  of  the  beam  and  surface-bonded  piezoceramic. 

The  ( t )  notation  has  been  dropped  for  convenience.  For  this  work  we  will  be  concerned 
with  the  steady-state  response  of  the  absorber,  therefore  we  transform  equation  (9)  into  the 
frequency  domain: 


El 


^  {Mb  +  fiMp) 


1  A  B 


(10) 


Actuator  Model 


The  actuators  are  modeled  with  the  one-dimensional  constitutive  relationships  for  a  piezo¬ 
ceramic  material: 


Dz  —  e33^3  +  d'32T,2 


e  —  d^Ez  +  I22T2 


(11) 


where  £33,  <^32 ,  and  Y22  are  material  properties  of  the  ceramic.  D3  is  the  electric  displacement, 
Ez  is  the  electric  field  applied  in  the  z  direction,  and  e  and  T2  are  the  strain  and  stress, 
respectively.  The  applied  electric  field  Ez  is  simply  the  applied  voltage  divided  by  the 
thickness  of  the  piezoceramic,  ha.  The  applied  stress  is  equivalent  to  the  actuator  force 
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> 


divided  by  the  side  area  of  the  ceramic,  Integrating  the  expression  for  the  electric 

displacement  over  x  and  y  yields  an  expression  for  the  induced  charge,  qa,  and  displacement, 
xa,  of  the  piezoceramic: 


Q  =  CV  +  x0fa 

xa  =  xcM  +  l rfa- 


(12) 


The  moment  induced  by  the  piezoceramic  actuators  can  be  written  as 


ma  =  2  fa 


hb  +  ha 
2 


xo^a\^  kaxa]  (^6  T  ^a) 


(13) 


Assuming  small  deformations  and  hence  small  rotation  angles,  we  can  write  the  extension 
of  the  actuator  as 


xa  = 


hb  +  ha 
2  L 


(14) 


Equations  (13)  and  (14)  can  be  combined  to  yield  an  expression  for  the  moment  applied 
by  the  piezoceramic  actuator  pair.  Transforming  the  expression  into  the  frequency  domain 
yields 

MaW  =  X0ka  (hb  +  ha)  V(u)  -  ka^^-BjQ  (u)  (15) 


Acoustic  Impedance  Model 

The  final  component  of  the  absorber  model  is  an  expression  for  the  acoustic  load  applied  to 
the  piston.  We  will  utilize  frequency-dependent  impedance  models  as  a  means  of  related  the 
piston  motion  to  the  sound  pressure.  The  use  of  impedance  models  makes  it  straightforward 
to  vary  the  type  of  acoustic  load  presented  to  the  piston.  A  general  impedance  model  of  the 
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acoustic  load  is 

W)=ZM'  (16) 

where  U  is  the  volume  velocity  of  the  piston.  Assuming  that  the  volume  velocity  is  equal 
to  Apw(l,t)  due  to  the  fact  that  the  displacement  of  the  piston  is  equivalent  to  the  tip 
displacement  of  the  beam,  we  can  write 


ApP(u)  =  juApZa  (u)  BpQ(uj) . 


(17) 


For  this  work  we  will  assume  that  the  absorber  is  connected  to  a  cylindrical  tube 
with  a  rigid  end  termination.  Under  this  assumption,  the  acoustic  impedance  is  Za(u>)  — 
—j  ( Pac/Ap )  cot(ul/c),  where  pa  is  the  density  of  air,  c  is  the  speed  of  sound,  and  l  is  the 
length  of  the  tube.  This  function  can  be  approximated  in  the  frequency  domain  as 


_  1  A^niUl-o;2/^) 

Arl  n" ,  (l  - 


1  Po<? ,  ,  , 

- - —kt[u). 

juApl  tK  1 


For  a  cylindrical  tube  with  a  rigid  termination: 


(18) 


uzi  =  (2 i  -  1)  f 

Wpj  =  if 


i  =  1,2,... 
^  =  1, 2 — 


Combining  equations  (17)  and  (19)  yields  an  expression  for  the  force  applied  to  the  piston 
as  a  result  of  flexure  motion: 


ApPM  =  (20) 

The  term  kt(ui)  is  a  frequency  dependent  stiffness  that  will  vary  depending  on  the  type  of 
acoustic  load  that  is  presented  to  the  absorber. 
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Table  1:  Nondimensional  parameters  for  the  absorber  analysis. 


Symbol 

Representation 

mass  ratio 

A* 

rrtpj  pAL 

actuator-to-flexure  stiffness  ratio 

AC  a 

e  Ea  L  (l+hb/haf 

E  La  (hb/haf 

tube-to-flexure  stiffness  ratio 

Kt 

Pa.C^  Ap /l 

EI/L 3 

structure-to-acoustic  time  scale 

7 

V  EI/pAL * 

nN  ( 1  .  4  n-2) 

tube  stiffness 

kt  (a) 

lli=l\^  72tt2(2i  — 1)^  J 

Nondimensional  Analysis 


Equations  (10),  (15),  and  (20)  are  combined  into  the  transducer  equations  of  motion  and 
nondimensionalized  for  the  parametric  analysis.  The  equations  are  nondimensionalized  by 
introducing  the  following  change  of  variables: 


o  =  t o^JpAlJ/EI  rj  =  q/L 


(21) 


The  resulting  nondimensional  equations  are 


Kb  +  KaBaBl  +  Ktkt{a)BpBj  -  a2  (Mb  +  pMp)\  p  (a)  =  2 Ka^~ - j——BaV  (a) .  (22) 

J  fla  1  +  /ib/fl'd 

The  terms  Ka,  Kt,  and  kt  (a)  are  defined  in  Table  1.  Defining 


H{a)  =  Kb  +  KaBaBTa  +  Ktkt(a)BpBp  -  a2  (Mb  +  pMp) 


(23) 


we  can  write 


x0v(ff)  ~  ko.^H  (a)  Ba 


,V(a) 

^p(g) 

x0V{a )  2  L 


Mulhsddn  =  KaBrH  Ba 


(24) 
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Figure  3:  (a)  Magnitude  of  ^~v(£)  f°r  three  values  of  fi:  1/10  (dotted),  1  (solid),  and  10 
(dashed);  (b)  magnitude  of  for  three  values  of  Ka:  1/1000  (dotted),  1  (solid),  1000 

(dashed). 

The  expressions  in  equation  (24)  represent  the  transfer  functions  between  the  applied  voltage 
and  the  actuator  motion,  xa,  and  the  tip  motion  of  the  flexure,  xp. 

Equation  22  illustrates  that  the  response  of  the  transducer  is  a  function  of  the  mass 
ratio,  the  actuator-to-flexure  stiffness,  the  tube  stiffness,  and  the  tube-to-flexure  stiffness 
ratio.  The  tube  stiffness  is  a  function  of  the  number  of  terms  in  the  acoustic  expansion 
and  the  structure-to-acoustic  time  scale.  The  effect  of  varying  these  four  parameters  is 
studied  with  a  numerical  model  of  the  transducer.  The  numerical  model  utilizes  a  five-mode 
expansion  of  a  clamped-sliding  beam  and  a  ten-mode  expansion  of  the  acoustic  impedance. 
The  mode  shapes  and  natural  frequencies  of  a  clamped-sliding  beam  can  be  found  in  Inman. 
For  brevity  they  will  not  be  repeated  here. 
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The  four  parameters  fi,  Ka,  Kt,  and  kt  (a)  were  studied  with  the  numerical  model  to 
understand  their  effect  on  the  transducer  response.  Figure  3a  is  a  plot  of  \Xa/x0V\  as 
a  function  of  the  mass  ratio  fj,.  As  expected,  increasing  the  mass  ratio  from  1/10  to  10 
increases  the  ratio  between  the  second  elastic  mode  of  structure  and  the  fundamental  natural 
frequency.  For  the  case  when  fi  =  1/10  the  frequency  ratio  is  approximately  5.2,  which  is 
the  frequency  ratio  a  bare  beam  with  clamped-sliding  boundary  conditions.  As  the  tip  mass 
is  increased  to  10  the  frequency  ratio  becomes  closer  to  20  and  the  flexure  begins  to  exhibit 
the  dynamics  of  a  single-mode  system. 

A  similar  analysis  is  performed  for  the  actuator-to-structure  stiffness  ratio.  When  the 
stiffness  ratio  is  much  less  than  1,  the  magnitude  of  the  actuator  motion  is  almost  four  orders 
of  magnitude  smaller  than  the  free  extension  of  the  actuator  (x0V).  Increasing  the  stiffness 
ratio  to  unity  increases  the  actuator  displacement  considerably  without  a  significant  effect 
on  the  pole-zero  spacing  in  the  frequency  response  \Xa/x0V\.  In  the  case  when  Ka  »  1,  the 
average  displacement  of  the  actuator  is  approximately  equal  to  the  free  displacement  x0V, 
but  the  mismatch  between  the  actuator  stiffness  and  the  flexure  stiffness  creates  a  pole-zero 
cancellation  in  the  frequency  response.  This  pole-zero  cancellation  is  due  to  the  negligible 
resistance  of  the  beam  to  actuator  motion. 

These  analyses  illustrate  variation  in  the  mass  ratio  and  actuator-to-structure  stiffness 
ratio  without  any  coupling  to  the  acoustic  field.  Structural-acoustic  coupling  is  introduced 
by  setting  nt  >  0  and  varying  the  tube-to-fiexure  stiffness  ratio.  Figure  4a  illustrates  the 
change  in  the  frequency  response  \Xa/x0V\  as  a  function  of  Kt.  When  Kt  «  1  there  is 
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(a)  (b) 


Figure  4:  (a)  Magnitude  of  for  two  values  of  Kt:  1/100  (dashed),  1  (solid)  . 

negligible  coupling  between  the  acoustic  field  and  the  actuator  motion.  Increasing  Kt  to  unity 
introduces  a  number  of  resonances  into  the  frequency  response.  These  resonances  correspond 
to  the  coupled  structural-acoustic  modes  of  the  enclosed  cavity.  Figure  4a  demonstrates  the 
basic  feasibility  of  coupling  the  piezoceramic  transducer  to  an  enclosed  cavity  modeled  as  a 
finite  tube.  The  coupling  between  the  acoustic  field  and  the  transducer  response  manifests 
itself  in  the  form  of  an  interlacing  pole-zero  pattern  of  structural-acoustic  resonances.  A 
beneficial  aspect  of  this  frequency  response  is  that  the  actuator  and  sensor  are  collocated 
with  one  another.  Sensor-actuator  collocation  enables  a  simple  and  robust  method  of  energy 
dissipation  in  the  resonances  through  feedback  control. 
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Feedback  Control  Simulations 


A  simple  method  of  estimating  the  effectiveness  of  feedback  control  is  to  simulate  active 
damping  as  a  means  of  energy  dissipation.  The  results  of  the  previous  section  illustrated 
that  structural-acoustic  coupling  is  possible  through  judicious  choice  of  transducer  design 
parameters.  Namely,  minimizing  the  mass  ratio  and  making  the  stiffness  ratios  approxi¬ 
mately  equal  to  unity  resulted  in  coupling  between  the  collocated  sensor-actuator  response 
of  the  piezoceramic  transducer.  These  results  will  be  the  basis  for  a  preliminary  analysis  of 
active  damping  of  the  coupled  modes. 

Feedback  control  simulations  are  performed  by  substituting 

V  (o)  =  Vd  (o)  ~  jo-^—Xa  {o )  (25) 

X0f^a 

into  equation  (22),  where  Vd{o)  is  a  disturbance  voltage  and  kv  is  a  scalar  feedback  gain. 
Combining  equations  (14)  and  (22)  yields  an  expression  for  the  closed-loop  matrix 

Hd  (o)  =  Kb  +  KaBaBTa  +  Ktkt(a)BpBp  +  jak :vBaBTa  -  a 2  {Mb  +  fxMp) ,  (26) 

which  can  be  directly  substituted  into  the  transfer  functions  in  equation  (24)  to  simulate  the 
closed-loop  response. 

Figure  4b  illustrates  the  effect  of  velocity  feedback  on  the  coupled  structural-acoustic 
modes.  The  design  parameters  chosen  for  the  simulation  are  /z  =  1/10,  Ka  =  Kt  =  l.  The 
simulation  of  the  closed-loop  response  demonstrates  the  broadband  damping  that  is  achieved 
with  velocity  feedback  of  the  piezoceramic  sensor  output.  The  first  eight  structural-acoustic 
modes  have  approximately  the  same  level  of  damping. 
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Summary  and  Future  Work 

The  critical  parameters  associated  with  utilizing  piezoceramic  sensor-actuators  for  acoustic 
control  were  derived  in  this  paper.  The  model  consisted  of  an  acoustic  source  driving  a 
cylindrical  cavity.  Four  nondimensional  parameters  were  identified:  the  mass  ratio  between 
the  acoustic  piston  and  the  beam  flexure,  the  actuator-to-flexure  stiffness  ratio,  the  acoustic- 
to-flexure  stiffness  ratio,  and  the  frequency  response  of  the  acoustic  load.  Numerical  analysis 
demonstrated  that  minimizing  the  mass  ratio  and  making  the  stiffness  ratios  approximately 
unity  produced  a  transducer  that  was  coupled  to  the  acoustic  field.  Feedback  control  sim¬ 
ulations  illustrated  the  ability  of  the  transducer  to  introduce  broadband  damping  into  the 
coupled  structural-acoustic  modes. 

Future  work  will  concentrate  on  quantifying  the  control  authority  of  the  transducer  and 
experimentally-verifying  transducer  performance.  A  quantitative  measure  of  the  control 
authority  will  be  studied  by  developing  a  state-space  model  of  the  transducer  coupled  to  the 
enclosed  cavity.  This  will  allow  us  to  use  measures  such  as  controllability  and  observability 
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grammians  to  quantify  the  coupling  of  the  transducer  to  select  modes.  It  will  also  allow  us 
to  perform  a  more  detailed  search  of  the  entire  parameter  space  and  introduce  optimization 
procedures  for  maximizing  the  structural-acoustic  coupling.  Once  these  studies  have  been 
performed  we  will  produce  a  prototype  transducer  and  experimentally  verify  the  control 
performance. 
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INVESTIGATION  INTO  THE  POWER  LOSSES  FROM  AMTEC  COMPONENTS 
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Abstract 

In  the  recent  past  a  number  of  programs  have  focused  on  developing  AMTEC  technology  associated  with  its 
fabrication  and  design.  The  performance  level,  however,  achieved  to-date  is  below  the  theoretical  potential 
of  this  device.  For  improving  performance  characteristics  we  examined  the  sources  of  power  loss  and 
proposed  some  solution  which  are  demonstrated  with  computer  calculations.  Two  kind  of  losses  that  reduce 
efficiency  are  thermal  (as  a  result  of  thermal  conduction  and  radiation  of  materials)  and  electrical  (related  to 
the  ohmic  resistance  of  material).  Each  of  the  losses  can,  in  principle,  be  reduced  separately  by  varying 
current  in  the  load.  To  reduce  the  thermal  losses,  the  current  must  be  increased;  to  reduce  the  electrical 
losses,  the  current  in  the  load  must  be  decreased.  In  such  inversely  competitive  situation  an  optimum  value 
must  be  sought  which  has  been  achieved  by  applying  the  optimization  theory.  Heat  losses  due  to  radiation 
are  reduced  by  increasing  the  current  density  and  by  reducing  the  emissivity  of  electrodes  and  other  surfaces. 
Changing  of  some  of  the  materials  used  in  the  cell  has  proved  to  be  helpful  in  improving  the  cell 
performance.  As  a  result  of  this  overall  effort  we  have  been  able  to  demonstrate  the  improvement  in  the 
efficiency  of  AMTEC  cell  by  77%. 
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Introduction 

A  power  system  considered  for  a  number  of  potential  space  missions  requires  long  life,  high  specific  power 
(power/mass),  that  is  low  mass,  high  areal  power  density  (power/  area),  high  efficiency,  low  cost  and  static 
in  operation.  To  a  great  extent  a  system  called  Alkali  Metal  Thermal  to  Electric  Converter  (AMTEC)  does 
possess  latent  qualities  that  would  qualify  it  to  be  a  possible  candidate  for  space  power.  It  can  provide 
efficiency  close  to  the  theoretical  Carnot  efficiency  at  relatively  low  temperatures.  AMTEC  has  conversion 
efficiency  much  higher  than  other  direct  thermoelectric  devices.  An  optimized  AMTEC  can  potentially 
provide  a  theoretical  efficiency  when  operated  between  1,000  K  to  1,300  K  on  the  hot  side  and  between  400 
K  to  700  K  on  the  condenser  side.  It  is  fuel  source  insensitive.  It  can  utilize  heat  as  input  fuel  from  most  of 
any  source  like  fossil  fuel,  the  Sun,  radioisotopes,  or  the  nuclear  reactor.  AMTEC,  with  solar  energy  as  a 
heat  source,  is  capable  of  being  an  alternative  to  photovoltaic-based  power  system  for  use  of  low  earth  orbit 
(LEO)  for  future  NASA  and  Air  Force  missions.  It  is  intended  to  be  used  for  future  NASA  missions  like 
Pluto  Express  (PX)  and  Europa  early  in  the  millennium  2000  with  radioisotope  decay  as  its  heat  input. 


Since  middle  sixties  a  number  of  programs  in  developing  the  operating  principle,  design  and  technology  of 
AMTEC  has  been  evolving  rapidly.  Kummer  and  Weber  demonstrated  the  conversion  of  heat,  through  the 
sodium  cycle,  into  electricity  by  the  use  of  beta" -alumina  solid  electrolyte  (BASE)  in  a  patent  assigned  to 
Ford  Motor  Company  in  1968  [1],  Several  years  later  Weber  described  the  operating  principle  of  AMTEC 
with  a  liquid  anode  in  a  historical  paper  in  1974  [2],  It  took  some  time  for  the  community  to  recognize  the 
potential  qualities  of  AMTEC  with  regard  to  its  uses  and  application  for  terrestrial  and  space  power. 
Accordingly,  there  was  not  much  activity  in  modeling,  design  and  technology  development  of  AMTEC  until 
almost  the  beginning  of  1990s.  During  that  AMTEC  dormant  period  there  were  few  landmark  papers 
primarily  on  its  principle  and  working  efficiency  [3-7].  A  review  article  on  the  progress  of  thermionic 
technology  during  a  ten-year  period  of  1983-92  covers  the  slow  moving  development  in  AMTEC  technology 
during  that  period  [8].  The  use  of  nuclear  power,  from  the  radioisotopes  or  nuclear  reactors,  for  the  space 
pointing  out  the  potential  use  of  AMTEC  perhaps  provided  another  motivation  in  looking  for  its  application 
[9-11].  A  flurry  of  AMTEC  activities  began  with  a  great  vigor  around  the  world  at  the  beginning  of  1990s. 
NASA  Lewis  Research  Center  (LeRC)  has  an  interest  in  developing  thermal  energy  storage  (TES)  for  the 
solar  dynamic  ground  test  demonstration  (SDGTD)  program.  The  research  and  development  effort  at  Jet 
Propulsion  Lab  (JPL)  included  studies,  which  address  both  overall  device  construction  and  investigation  into 
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the  AMTEC  components  [12-25].  Advanced  Modular  Power  Systems  (AMPS)  have  been  focusing  on 
designing  and  manufacturing  AMTEC  cells  for  Pluto  Express  mission  [26-27].  Air  Force  Research  Lab 
(AFRL)  has  an  interest  in  developing  an  electrical  power  system  for  the  payload  in  low  earth  orbit  (LEO) 
for  a  duration  of  five  years  and  in  testing  the  performance  of  AMTEC  for  PLUTO  express  [28-30].  Orbital 
Science  Corporation  has  been  conducting  a  series  of  studies  of  radioisotopes  power  system  based  on  general- 
purpose  heat  source  (GPHS)  for  potential  deep  space  missions  [9,31-33].  AMTEC  has  been  found 
compatible  with  direct  conversion  from  the  projected  SP-10  nuclear  space  power  reactor  [10-11],  Recently 
University  of  New  Mexico  has  been  engaged  in  modeling  and  analysis  of  AMTEC  performance  and 
evaluation  [34],  In  Japan  the  interest  in  AMTEC  analysis  has  evolved  at  Kyushu  University  in  collaboration 
with  Electrotechnical  Laboratory  [35-36].  Creare  Inc.  has  primarily  concentrated  on  evaporator  component 
of  AMTEC  to  enable  it  operate  efficiently  [37-39]. 


While  these  programs  resolved  a  number  of  key  technological  issues  associated  with  the  design  and 
fabrication  of  AMTEC  successfully,  the  performance  level  achieved  hitherto  is  still  below  the  theoretical 
potential  of  this  device.  In  tracking  down  the  source  of  power  loss  during  the  test  cell  performance  of 
various  AMTEC  designs  the  major  losses  are  found  to  be  radiative  and  conductive.  Since  the  cell  is  tested  in 
vacuum  there  is  no  convective  loss.  Most  of  the  heat  escapes  through  the  cell  wall.  There  are  radiative 
losses  due  to  condenser  and  along  the  artery  as  well.  In  order  to  minimize  these  losses  we  propose  to  use 
some  different  materials  which  should  be  compatible  with  the  operating  conditions  of  the  cell  and  capable  of 
reducing  these  losses.  We  select  a  test  AMTEC  cell  designated  by  PX-3A,  which  has  been  in  operation  at 
US  Air  Force  Research  Laboratory  since  July  9,  1997  and  compare  its  theoretical  performance  as  a  function 
of  the  variation  in  the  material  parameters. 


Principle  and  Working  of  AMTEC  Cell 

A  typical  AMTEC  cell  is  shown  in  Fig.  1  schematically.  An  AMTEC  cell  is  a  static  device  for  the 
conversion  of  heat  to  electricity  performing  two  distinct  cycles:  1)  conversion  of  heat  to  mechanical  energy 
via  a  sodium-based  (or  any  suitable  alkali  metal)  heat  engine  and  2)  conversion  of  mechanical  energy  to 
electrical  energy  by  utilizing  the  special  properties  of  the  beta"  aluminum  solid  electrolyte  (BASE)  material. 
It  has  the  potential  for  reliable  long  life  operation  and  high  conversion  efficiency.  The  general  principles 
governing  the  operation  of  an  AMTEC  cell  have  been  given  quite  elaborately  in  the  early  stages  of  its 
development  [1-5].  The  essential  thermodynamic  aspects  of  the  sodium-based  engine  have  been  described  in 
those  references.  Accordingly  a  brief  outline  will  be  presented  here  for  the  self-consistency  of  this  work. 
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AMTEC,  a  relatively  new  type  of  device,  is  based  on  the  principle  of  a  sodium  concentration  cell,  conceived 
in  late  sixties  [1-4].  A  closed  vessel  is  divided  into  a  high  temperature,  high-pressure  region  in  contact  with 
a  heat  source  and  a  low  temperature,  low  pressure  region  in  contact  with  a  heat  sink.  A  barrier  of  a  BASE 
sheet  whose  ionic  conductivity  is  much  larger  than  its  electronic  conductivity  separates  these  regions.  The 
BASE  is  coated  with  a  porous  metal  electrode  (cathode)  which  covers  the  low-pressure  surface  of  the  BASE. 
A  closed  container  is  partially  filled  with  a  small  quantity  (typically  <10  g)  of  liquid  sodium  as  the  working 
fluid.  Sodium  ions  disperse  through  the  BASE  in  response  to  the  pressure  differential  (gradient  of  Gibbs  free 
energy).  Electrical  leads  are  connected  with  the  porous  cathode  and  with  the  high  temperature  liquid  (vapor) 
sodium,  which  acts  as  the  anode.  When  the  circuit  is  closed  electrons  flow  to  the  porous  anode  surface 
through  the  load,  producing  electrical  work.  A  return  line  and  electromagnetic  pump  or  a  wick  circulates  the 
sodium  from  the  cold  zone  (condenser)  to  the  hot  zone  (heater)  of  AMTEC.  It  is  limited  to  1350  K  because 
of  sodium  interaction  with  the  BASE.  The  lower  temperature  region  is  limited  to  a  minimum  of  500  K  by 
the  need  of  maintaining  the  sodium  in  the  liquid  state  (but  it  would  be  operating  around  623  K  for  high 
efficiency).  As  the  condenser  temperature  increases  above  700  K  the  efficiency  decreases  [3]. 

The  main  component  of  an  AMTEC  cell  is  the  BASE  tube  which  essentially  conducts  sodium  ions  much 
more  rapidly  than  what  it  does  for  neutral  sodium  atoms  or  negative  particles  like  electrons.  A  sodium 
pressure  difference  across  a  thin  BASE  tube  drives  sodium  ions  from  the  high-pressure  side  to  the  low- 
pressure  side.  Thus  positively  charged  sodium  ions  concentrate  on  the  low-pressure  side  while  negatively 
charged  particles,  electrons  remain  concentrated  on  the  high-pressure  side  (inside  the  BASE  tube),  resulting 
an  electrical  field  gradient  across  the  BASE  tube.  This  gradient  balances  the  sodium  pressure  differential 
thus  preventing  the  further  flow  of  sodium  ions.  This  electrical  potential  difference  is  utilized  to  drive  an 
electric  current  through  a  load  by  applying  some  appropriate  electrodes  at  the  two  ends  of  the  tube.  The 
unusual  properties  of  the  BASE  provide  a  method  of  converting  the  mechanical  energy,  represented  by  the 
pressure  differential,  into  electrical  energy.  Or  in  other  words  the  BASE  converts  a  chemical  potential 
difference  into  an  electrical  potential  difference.  As  the  electrons,  after  going  through  the  load,  reunite  with 
sodium  ions  they  form  neutral  sodium  atoms  in  the  vapor  state  on  the  low-pressure  side.  The  sodium  vapor 
travels  to  the  condenser  where  it  condenses  into  liquid  state.  The  sodium  liquid  is  pressurized  through  a 
wick  or  electromagnetic  pump  and  is  returned  to  the  high-pressure  side  of  the  BASE.  In  that  way  the 
thermal-to-mechanical-to-electrical  conversion  process  is  completed.  The  efficiency  of  this  final  conversion 
is  governed  by  variety  of  irreversible  kinetic  and  transport  processes  occurring  at  the  electrode  interfaces, 
within  the  BASE  material,  internal  impedance,  and  thermal  conduction  and  radiation  losses  [12-14],  A 
number  of  efforts  are  underway  to  develop  practical  and  high  efficiency  cells  [18,  40-45].  AMPS,  in 
collaboration  with  AFRL,  has  been  manufacturing  and  testing  for  improving  efficiency  of  PX  series  of 
AMTEC  cells  [26-30].  The  latest  multi-tube  AMTEC  cells  have  been  in  operation  for  some  time.  The 
longest  in  operation  among  them  is  PX-3A  whose  parameters  are  given  in  Table  1.  We  have  modified  it 
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theoretically  for  improving  efficiency  by  varying  parameters  governing  its  material  properties  for  optimum 
conditions  in  a  computer  model  of  the  cell  to  predict  its  performance. 


AMTEC  Theory  and  Analysis 

The  transport  of  sodium  through  an  AMTEC  cell  is  a  complex  phenomenological  process.  Its  exact  detailed 
analysis  is  very  difficult  and  is  further  complicated  as  it  requires  the  simultaneous  solution  of  thermal,  fluid 
flow,  and  electrical  equations.  Those  equations  are  interdependent,  with  each  of  the  three  analyses  requiring 
the  results  of  the  other  two.  This  interdependence  is  between  a  number  of  axially  varying  distribution 
functions.  (Had  it  been  between  single- valued  variables ;  it  could  have  been  solved  by  a  set  of  simultaneous 
algebraic  equations.)  Specifically  ,  solving  for  the  cell's  temperature  distribution  requires  the  knowledge  of 
the  axial  variation  of  the  sodium  flux  through  the  BASE  tube  and  of  the  electrical  output  power  density 
profile  over  the  tube  length.  Solution  of  the  axial  pressure  variation  of  the  low-pressure  sodium  requires 
knowledge  of  the  cell’s  temperature  distribution  and  the  BASE  tubes’  current  density  variation.  Similarly, 
solving  for  the  axial  variation  of  the  current  density  and  of  the  inter-electrode  voltage  requires  prior  solution 
of  the  axial  variation  of  the  BASE  tube  temperature  and  intemal-to-extemal  pressure  ratio.  Those 
interdependent  distribution  functions  require  solutions  of  coupled  differential  and  integral  equations  by  a 
more  sophisticated  procedure.  Schock  et  al  [32]  generated  a  thermal  analysis  model  for  multi-tube  AMTEC 
cell  by  appropriately  modifying  the  IT AS  and  SINDA  codes  [46,47], 


Energy  conversion  devices  have  few  equilibria  and  are  typically  open  systems  unlike  the  classical 
thermodynamics  which  is  restricted  to  reversible  and  closed  systems.  Onsager’s  treatment  of  irreversible 
processes,  such  as  diffusion,  can  be  applied  to  AMTEC  operation  to  deal  with  its  irreversibility  of  the 
process  and  openness  of  the  system.  One  can  write,  in  principle,  the  effective  emf,  V  as  a  function  of  cell 
voltage  in  open-circuit,  Voc  and  electrode  polarization  over  potential,  V0p  =  (^a  -  )  [2]  from  Nemst 

equation  [5].  The  open-circuit  voltage  and  charge-exchange  current  density  are  related  with  the  cell 
temperature  and  pressure  [19].  We  can  thus  write  for  the  net  power  output  as: 


Pout  -  VI, 


(i) 


where 


V  -  Vqc  -  JRint  -  V0p,  (2) 

Voc  =  RTfi/F  In  Pa/Pc>  (3) 


R  is  the  gas  constant  =  8.314  J/mol.K.  Tg  is  the  temperature  of  the  BASE  tube,  Pa  and  Pc  are  the  pressures 
at  the  anode  and  cathode  sides  respectively.  I  is  the  net  current  in  the  circuit  and  Rint  is  the  total  internal 
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resistance  of  the  cell  including  the  contact  resistance  of  the  electrodes  Rcontact>  sheet  resistance  in  the  plane 
of  electrodes  Rsheet.  resistance  of  the  current  collectors  Rcollector.  rsistance  of  the  bus  wires  and  conductor 
leads  to  the  load  Rbus>  the  resistance  for  charge  exchange  polarization  losses  at  the  BASE-electrodes 
interfaces  Rop,  and  the  BASE  ionic  resistance,  Rb  given  by: 

RB=PBtB  (4) 


where  p£and  tB  are  electrical  resistivity,  the  expression  for  which  has  been  developed  by  Steinbruck  [49], 

and  the  thickness  of  BASE  tube  respectively.  Rjnt  will  be  written  here  as  the  sum  of  the  terms  in  the  same 
order  as  stated  above: 


Heat  Losses 


Rint  =  Rcontact  +  Rsheet  +  Rcollector  +  Rbus  +  Rop  +  Rb  (5) 


For  a  multitube  AMTEC  cell,  the  heat  losses  are  shown  in  Figure  1.  The  most  of  the  heat  passes  the  cooling 
fluid  (air)  through  the  condenser.  The  heat  passes  the  cooling  fluid,  Qair,  is  result  of  heat  conduction,  heat 
radiation  and  condensation  of  vapor  sodium  as: 

Qair  =  Qwall  +  Qartery  +  Qrcond  +  Qcond  -  Qcold  (6) 

where  Qwall  is  the  conducted  heat  to  condenser  through  the  cell  wall,  Qartery  is  conducted  heat  through  the 
artery,  Qrcond  is  net  radiated  energy  between  the  condenser  and  other  surfaces,  Qcond  is  latent  heat  of 
condensation  of  sodium,  and  Qcold  is  heat  conduction  loss  at  edge  of  cold  plate.  These  heat  losses  can  be 
expressed  in  the  following  expressions: 

Qwall =  kwall  A-wall  (dT/dx)x=at  condenser  (7) 

Qartery  =  kartery  A  artery  (dT/dx)x=at  condenser  (8) 

Qcond  =  hfgWl  (9) 

where  k  is  thermal  conductivity,  A  is  area,  hfg  is  latent  heat  of  condensation  of  sodium  per  unit  mass  and  ril 
is  sodium  flow  rate.  For  wall  and  artery  dT/dx  is  the  temperature  gradient  at  the  condenser.  As  an  enclosure, 
every  surface  in  the  cell  effects  heat  radiation.  Net  radiated  energy  between  the  condenser  and  other 
surfaces  depends  on  view  factors  Fjj,  emissivities  ej,  and  temperatures  Tj  for  every  surfaces,  so  it  can  be 
expressed  as 

Qrcond  =  ffcl,  £2£3. . ,Tl,  T%  T3, . ,Fn,  Fi2,Fi3, . )  (10) 

Some  heat  goes  to  ambient  through  cell  wall  and  insulation;  and  it  can  be  given  as 
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x=cell  length 

Qinsul  —  J  K(x)P[Twan  (x)  —  TinSui  (x)]lx  (11) 

x=0 

where  K(x)  is  some  constant  which  includes  conduction  and  radiation  effects,  and  P  is  the  perimeter  of  cell 
wall.  Twan  and  Tjnsui  are  cell  wall  temperature  and  insulation  temperature  which  change  along  the  cell 
length.  Heat  conduction  loss  at  edge  of  cold  plate,  Qcold  .  can  be  calculated  similar  to  Qinsul  taking 
temperature  TCOnd  instead  of  TWall,  but  Qcold  is  small.  The  total  of  heat  losses  become 

Qloss =  Qair  +  Qinsul  +  Qcold  (12) 

The  overall  conversion  efficiency  of  AMTEC  cell  is  given  by  [5,50]: 

"H  =  Pout/Qinput  =  VI  /  [VI  +  Qloss]  (13) 

In  order  to  get  the  maximum  efficiency  the  total  heat  losses,  Qloss  .  must  be  minimum.  This  implies  that: 

•  The  choice  of  working  fluid  be  such  that  its  latent  heat  of  condensation  must  be  low 

•  The  emissivities  of  materials  used  for  BASE  tube,  outer  condenser  surface,  outer  artery  surface,  inner 
cell  wall  and  radiation  shields  must  be  low  for  minimizing  radiation  losses 

•  For  conductive  losses  through  the  wall,  artery  and  insulation,  low  thermal  conductive  materials  must  be 
used,  but  for  hot  plate  high  thermal  conductivity  material  must  be  used. 

•  The  output  power  maximizes  when  the  load  resistance  is  equal  to  the  internal  resistance  of  the  cell.  This 
implies  that  the  load  resistance  should  match  the  internal  resistance  for  the  best  results. 

We  then  proceeded  with  these  constraints  as  input  data  for  the  computer  program  in  order  to  get  the  resulting 
output  power  and  efficiency 

Results 

I  consider  the  parametric  properties  of  materials  in  our  optimization  analysis.  The  effects  of  material  on 
electrical  power  generated  by  the  PX-3A  cell,  the  amount  of  the  heat  supplied  to  the  cell,  the  efficiency  of 
the  cell,  and  the  losses  from  the  cell  are  given  in  figures  2-9.  The  temperature  maintained  at  the  hot  plate  is 
1173  K,  and  at  the  condenser  side  is  623  K  for  all  these  cases.  In  these  figures  the  dotted  line  shows  the 
values  for  original  PX-3A  cell.  The  maximum  power  generated  by  the  cell  is  given  in  Table  2.  Some  other 
relevant  parameters  for  the  maximum  power  are  given  in  the  same  table.  Also  heat  losses  are  given  as  the 
percentage  of  the  Qinput-  The  following  material  changes  are  made  for  the  cell  to  compare  the  results: 


Case  A:  material  for  hot  plate,  stud  and  evaporator  is  changed  from  SS  to  Ni 

Case  B:  in  addition  to  case  A,  the  material  for  condenser,  artery,  shield  and  cell  wall  is  changed  from  SS  to 
Inconel 

Case  C:  in  addition  to  case  B,  the  emissivity  of  sodium  film  is  changed  from  0.05  to  0.025;  because  the 
emissivity  of  liquid  sodium  is  given  as  0.02  in  [51] 

Case  D:  in  addition  to  case  C,  cell  wall  and  artery  coated  with  Rh  (e=0.06) 

Case  E:  in  addition  to  case  D,  material  for  cell  wall  above  the  BASE  tube  is  changed  from  SS  to  ceramic, 
Zr02 

Case  F:  in  addition  to  case  E,  perfect  insulator  is  assumed  as  a  insulation  material  instead  of  Molded  Mink- 
K. 


In  Fig.  2,  it  can  be  seen  that  the  electrical  power  generated  by  the  cell  increases  appreciably  for  case  A.  In 
this  case  the  maximum  power  reaches  from  4.49  W  to  5.97  W.  For  cases  B,  C,  D,  E  and  F,  the  electrical 
power  increases  slowly.  The  maximum  power  in  case  F  is  obtained  as  6.26  W.  Heat  input  to  the  cell  versus 
current  is  given  in  Fig.  3.  The  maximum  heat  goes  into  the  cell  in  case  A  because  Ni  has  higher  thermal 
conductivity  than  SS.  For  case  B,  the  amount  of  heat  supplied  to  the  cell,  Qinput.  decreases  because  of  the 
Inconel’s  lower  thermal  conductivity  which  impeds  the  loss  of  heat  through  walls  and  condenser.  Depending 
on  heat  losses,  Qinput  g°es  down  for  cases  C,  D,  E  and  F;  and  it  takes  the  minimum  values  in  case  F. 


The  cell  conversion  efficiency  versus  current  is  given  in  Fig.  4.  The  maximum  efficiency  is  12.98%  for 
original  cell,  the  maximum  efficiencies  for  cases  A  through  F  are  substantially  increased  around  the  current 
2-2.5  A.  They  are  given  in  Table  3.  Heat  losses  through  cell  wall  and  insulation,  Qinsub  through  artery, 
Qartery.  and  through  cell  wall,  Qwall.  arc  given  in  Figures  5-7.  They  show  higher  values  from  original  cell 
in  cases  A,  B  and  C,  but  smaller  values  in  cases  D,  E  and  F.  Heat  loss  in  case  F  is  minimum  for  Qinsul.  but 
in  case  of  Qartery  and  Qwall  the  minimum  is  for  case  E.  In  case  F,  the  perfect  insulator  is  assumed  for 
insulation  material.  It  means  there  is  no  heat  losses  through  cell  wall  and  insulation.  That  time  temperatures 

of  the  cell  elements  (wall,  evaporator,  BASE  tube,  artery, . )  become  higher  (Table  3).  Therefore  Qartery 

and  Qwall  will  increase  when  Qinsul  become  zero.  As  the  conducted  heat  losses  through  artery  for  cases  D 
and  E  are  almost  the  same,  the  conducted  heat  losses  through  the  cell  wall  for  cases  D  and  E  show  big 
difference.  Because  the  ceramic  material  Zr02  is  considered  only  for  cell  wall  not  for  artery,  in  case  E.  The 
Zr02  has  lower  thermal  conductivity  than  SS  has  (for  ZrC>2  2.5  W/mK  and  for  SS  25  W/mK  at  1000  K).  The 
net  heat  radiation  to  condenser  increases  for  cases  A  and  B,  but  it  decreases  for  other  cases  comparing  to 
original  cell.  Depends  on  the  low  sodium  film  emissivity  (e=0.025),  it  takes  the  smallest  values  in  case  C. 


Although  the  sodium  film  emissivity  is  the  same  for  cases  D,  E  and  F,  Qrcond  increases  because  of  higher 
temperatures  of  cell  elements  (Fig.  8).  The  condensation  heat  losses,  Qcond  are  almost  same  values  in  all 
cases  as  should  be  expected,  see  Fig.  9. 


Conclusion 

We  have  been  able  to  reduce  the  losses  almost  on  all  counts  thereby  increasing  the  output  power  and  thus  the 

efficiency  of  AMTEC.  The  efficiency  has  improved  considerably  when  we  incorporate  the  variation  in  the 

material  used  in  the  cell.  Besides  we  make  the  following  recommendations: 

•  The  interior  cell  wall  could  be  coated  with  some  appropriate  material  with  high  reflective  coefficient 
and  low  emissivity 

•  The  cell  wall  could  be  made  of  double  layers  with  air  in  between 

•  The  artery’s  outer  surface  should  be  coated  with  the  same  material  used  for  the  cell  wall  coating 

•  Condenser  may  be  coated  with  sodium  also 

•  The  power  output  maximizes  at  a  certain  value  of  the  current  in  the  circuit,  and  when  the  load 
resistance  matches  the  internal  resistance  of  cell.  For  best  results  the  cell  must  be  operated  under 
those  conditions. 

•  For  conductive  losses  through  the  wall,  artery  and  insulation,  low  thermal  conductive  materials  must  be 
used,  but  for  hot  plate  high  thermal  conductivity  material  is  recommended. 

•  A  PX-3A  AMTEC  cell  has  been  in  operation  at  US  Air  Force  Research  Lab,  Albuquerque  NM  since 
July  9,  1997  [52].  It  had  about  16%  degradation  in  performance  in  less  than  first  500  hours  of 
operation.  This  degradation  is  attributed  to  the  maturing  of  electrodes.  The  cell  degradation  continues 
after  the  electrode  maturing  period,  though  at  a  slower  pace,  see  Fig.  10.  The  cell  performance  can,  in 
principle,  be  enhanced  by  decreasing  the  parasitic  heat  losses. 


Acknowledgment 

I  am  indebted  to  Clay  Mayberry  and  John  Meril  for  Providing  the  AMTEC  data  before  publication  and  Dr. 
Jean-Michel  Toumier  for  many  useful  discussions.  This  work  is  based  in  parts  supported  by  AFOSR  Sub¬ 
contract  99-0832  CFDA  #12.800 


5-10 


References 


1.  J.T.Kummer  and  N.  Weber,  U.S.  Patent  3,458,356,  Assigned  to  Ford  Motor  Co,  (1968). 

2.  N.  Weber,  Energy  Conversion, 14, 1  (1974). 

3.  T.  K.  Hunt,  N.  Weber  and  T.  Cole,  in  Proceeding  of  the  13  th  Inter  society  Energy  Conversion 
Engineering  Conference,  SAE,  p.  2001,  Warrandale,  PA  (1978). 

4.  Hunt,  T.  K.,  N.  Weber,  and  T.  Cole,  in  Solid  State  Ionics,  J.  B.  Bates  and  G.  C.  Ferrington,  Editors,  p. 
263,  North  Holland  Pub.  Co.,  Amsterdam  (1981). 

5.  T.  Cole,  Science,  221, 915  (1983). 

6.  R.  Ewell  and  J.  Mondt,  in  Space  Nuclear  Power  Systems  1984,  M.  S.  El-Genk  and  M.  D.  Hoover, 
Editors,  p.  385,  Orbit  Book  Co.,  Malabar,  FL  (1985). 

7.  C.  P.  Bankston,  T.  Cole,  S.  K.  Khanna,  and  A.  P.  Thakoor,  in  Space  Nuclear  Power  Systems,  M.  S.  El- 
Genk  and  M.D.  Hoover,  Editors,  p.  393,  Orbit  Book  Co.,  Malabar,  FL  (1985). 

8.  Dahlberg,  R.  C.,  et  al,  "Review  of  Thermionic  Technology  1983  to  1992,"  "A  Critical  Rev  of  Space 
Nuclear  Power  and  Propulsion  1984  -1992",  Ed.  M.  S.  El-Genk  (AIP  Press,  NY  1994)  pp  121-161 

9.  A.  Shock,  in  Proc.  15th  Intersociety  Energy  Conversion  Engineering  Conference,  Vol.  2,  p.  1032, 

AIAA,  NY  (1980). 

10.  D.  Birden  and  F.  A.  Angelo,  in  Proc.  18th  Intersociety  Energy  Conversion  Engineering  Conference,  p. 
61,  AICHE,  N.Y.  (1983). 

11.  E.  J.  Brat  and  G.  O.  Fitzpatrick,  Direct  Conversion  Nuclear  Reactor  Space  Power  Sy.,  Report  No. 

AFW ALTR-82-2073 ,  Vol  1 

12.  R.  M.  Williams,  B.  Jeffries-Nakamura,  M.  Underwood,  B.  Wheeler,  M.  Loveland,  S.  Kikkert,  J.  Lamb, 
T.  Cole,  J.  Kummer,  and  C.  P.  Bankston,  J.  Electrochem.  Soc.,  136,  893  (1989). 

13.  R.  M.  Williams,  M.  E.  Loveland,  B.  Jeffries-Nakamura,  M.  L.  Underwood,  C.  P.  Bankton,  H.  Leduc  and 
J.  T.  Kummer,/.  Electrochem.  Soc.  137,  1709  (1990). 

14.  R.  M.  Williams,  B.  Jeffries-Nakamura,  M.  L.  Underwood,  C.  P.  Bankton  and  J.  T.  Kummer,  J. 
Electrochem.  Soc.  137,  1716  (1990). 

15.  M.  A.  Ryan,  B.  Jeffries-Nakamura,  R.  M.  Williams,  M.  L.  Underwood,  D.  O'Connor,  and  S.  Kikkert,  in 
Proc.  26th  Intersociety  Energy  Conversion  Engineering  Conference,  American  Nuclear  Society,  5,  463 
(1991). 

16.  M.  A.  Ryan,  R.  M.  Williams,  B.  Jeffries-Nakamura,  M.  L.  Underwood,  and  D.  O'Connor,  JPL  New 
Technology  Report ,  18620-8166  (1991). 

17.  M.  A.  Ryan,  B.  Jeffries-Nakamura,  D.  O’Connor,  M.  L.  Underwood,  and  R.  M.  Williams,  in  Proc. 
Symposium  on  High  Temperature  Electrode  Materials  and  Characterization,  D.  D.  Macdonald  and  A.  C. 
Khanker,  Editors,  Electrochem.  Soc.  91-6, 115  (1991). 


5-1] 


18.  C.  B.  Vinning,  R.  M.  Williams,  M.  L.  Underwood,  M.A.  Ryan,  and  J.W.  Suiter,  in  Proc.  27th 
Intersociety  Energy  Conversion  Engineering  Conf.,  Society  of  Automobile  Engineers,  3 , 123,  (1992). 

19.  M.  L.  Underwood,  R.  M.  Williams,  M.  A.  Ryan,  B.  Jeffries-Nakamura,  and  D.  O'Connor,  in  Proc.  9th 
Symp.  On  Space  Nuclear  Power  Systems,  M.  S.  El-Ghenk  and  M.  D.  Hoover,  Editors,  American  Inst. 
Phys.3,  1331  (1992). 

20.  M.  L.  Underwood,  D.  O’Connor,  R.  M.  Williams,  B.  Jeffries-Nakamura  and  M.  A.  Ryan,  in  Proc.  27th 
IECEC,  3, 197  (1992). 

21.  M.  L.  Underwood,  B.  Jeffries-Nakamura,  D.  O'Connor,  M.  A.  Ryan,  J.  W.  Suitar,  and  R.  M.  Williams, 
in  Proc.  28th  IECEC,  1,  855  (1993). 

22.  M.  A.  Ryan,  R.  M.  Williams,  B.  Jeffries-Nakamura,  M.  L.  Underwood,  and  D.  O'Connor,  NASA  Tech. 
Briefs,  p.  80,  (1993). 

23.  M.  A.  Ryan,  R.  M.  Williams,  C.  Spaipetch,  A.  Kisor,  D.  O'Connor,  M.  L.  Underwood,  and  B.  Jeffries- 
Nakamura,  CONF940101,  AIP  Press,  NY,  p.  1495  (1994). 

24.  M.  A.  Ryan,  A.  Kisor,  R.  M.  Williams,  B.  Jeffries-Nakamura,  and  D.O'Connor,  in  Proc.  29th  IECEC,  2, 
877  (1994). 

25.  M.  A.  Ryan,  R.  M.  Williams,  M.  L.  Phillips,  L.  Lora,  and  J.  Miller,  in  Proc.  Space  Techology  and 
Applications  International  Forum-1998,  M.  S.  El-Ghenk,  Editor,  p.1607  (AIP  Press,  NY),  (1998). 

26.  R.  K  Sievers,  T.  K.  Hunt,  J.  F.  Ivanenok,  J.  Pantolin,  and  D.  A.  Butkiewicz,  in  Proc.  !0th  Space  Nuclear 
Power  and  Propulsion,  CONF-930103,  M.  S.  El-Genk  and  M.  D.  Hoover,  Editors,  p.  319  (AIP  Press, 
NY)  (1993). 

27.  R.  K.  Sievers,  J.  R.  Rasmussen,  C.  A.  Barkowski,  T.  J.  Hendricks,  and  J.  E.  Pantolin,  in  Proc.  STAIF-98, 
CONF-980103,  p.1479  (AIP  Press,  NY),  (1998). 

28.  M.  J.  Schuller,  P.  Haugen,  E.  Reiners,  J.  Merrill,  R.  K.  Sievers,  R.  Svedberg,  J.  F.  Ivanenok,  III,  C.  J. 
Crowley,  and  M.  G.  Izenson,  in  Proc.  31st  IECEC,  IEEE,  paper  no.  96184,  2,  877  (1996). 

29.  J.  Merrill,  M.  J.  Schuller,  R.  K.  Sievers,  C.  A.  Borkowski,  L.  Huang,  and  M.  S.  El-Genk,  in  Proc.  32nd 
IECEC,  IEEE,  paper  no.  97379,  2,  1184  (1997). 

30.  J.  M.  Merrill,  M.  Schuller,  and  L.  Huang,  in  Proc.  STAIF,  M.  S.  El-Genk,  Editor,  CONF-980103,  (AIP 
Press  NY)  1613,  (1998). 

31.  A.  Schock,  in  Proc.  15th  Intersoc.  Energy  Conversion  Engineering  Conf.,  1032  (1980). 

32.  A.  Schock,  H.  Norivan,  C.  Or,  and  V.Kumar,  in  Proc.  48th  Inti.  Astronomical  Congress,  1  (1997) 


5-12 


33.  A.  Schock,  H.  Norivan,  C.  Or,  and  V.Kumar,  Preprint  for  presentation  in  33fd  Intesoc.  Energy 
Conversion  Engineering  Conf.  (1998) 

34.  J.  M.  Toumier  and  M.  S.  El-Ghenk,  in  Proc.  15th  STAIF,  (AIP  Press  NY)  1576  (1998). 

35.  A.  Kato,  Nakata,  and  K.  Tsuchida,  in  Proc.  27th  IECEC,  American  Chemical  Soc.,  Paper  No.  929006, 
p.  1,(1992). 

36.  A.  Kato,  Nakata,  and  K.  Tsuchida,  in  Proc.  28th  IECEC,  American  Chemical  Soc.,  Paper  No.  93027,  p. 
809,  (1993). 

37.  C.  J.  Crowley  and  M.  G.  Izenson,  in  Proc.  10th  Symp.  On  Space  Nuclear  Power  and  Propulsion,  CONF- 
930103,  M.  S.  El-Genk,  Editor,  (AIP  Press,  NY),  897  (1993). 

38.  M.  G.  Izenson  and  C.  J.  Crowley,  in  Proc.  28th  Intersociety  Energy  Conversion  Engineering 
Conference,  Paper  no.  93221,  American  Chemical  Society,  1:  p.  829,  (1993). 

39.  M.  G.  Izenson  and  C.  J.  Crowley,  in  Proc.  31st  Intersociety  Energy  Conversion  Engineering 
Conference,  Paper  no.  96262,  IEEE,  4,  p.  2226  (1996). 

40.  M.  Sayer,  M.  F.  Bell,  and  B.  A.  Judd,  J.  Applied  Physics,  67,  832  (1990) 

41.  M.  L.  Underwood,  D.  O’Connor,  R.  M.  Williams,  M.  A.  Ryan,  and  C.  P.  Bankston,  /.  Propulsion  and 
Power,  8,  878  (1993). 

42.  M.  L.  Underwood,  R.  M.  Williams,  M.  A.  Ryan,  and  B.  Jeffries-Nakamura,  in  Proc.  Space  Nuclear 
Power  and  Propulsion,  CONF-  930103,  M.  S.  El-Genk  and  M.  D.  Hoover,  Editors,  p.  885,  AIP,  (1993). 

43.  A.  Schock,  H.  Noravian,  V.  Kumar,  and  C.  Or,  in  Proc.  32nd  Interdisciplinary  Energy  Conversion 
Conf,  #97529  ,2,1136  (1997). 

44.  M.  A.  K.  Lodhi,  Michael  Schuller,  and  Paul  Housgen,  in  Proc.  Space  Techology  and  Applications 
International  Forum,  M.  S.  El-Genk,  Editor,  p.  1285,  (AIP  Press,  NY),  (1996). 

45.  M.  S.  El-Ghenk,  and  J-M.  Toumier,  in  Proc.  5th  ESPC-98, 416  ,  p.  257,  (1998). 

46.  H.  Noravian,  in  Proc.  26th  International  Conference  on  Environmental  Systems,  No.  961376,  (1996). 

47.  J.  Gaski,  S INDA  (System  Improved  Numerical  Differencing  Analyzer),  vers.  1.315  from  Network 
Analysis  Associate,  Fountain  Valley,  CA,  (1987). 

48.  A.  M.  Straus,  and  S.  W.  Peterson,  in  Proc.  STAIF,  p.1571,  (1998). 

49.  M.  Steinbruck,  V.  Heinzel,  F.  Huber,  and  W.  Pepler,  in  Proc.  28th  IECEC  1. 1, 799,  (1993). 

50.  C.  J.  Crowly,  M.G.  Izenson,  P.  N.  Wallis,  R.  K.  Sievers,  and  J.  F.  Ivanenok  III,  in  Proc.  29th  IECEC, 
Vol.l,  p.  882,  AIAA,  (1994). 

51.  J.  M.  Toumier,  M.  S.  El-Ghenk,  M.  Schuller,  and  P.  Hausgen,  in  Proc.  14th  STAIF-97, 397,  (1997). 

52.  J.  M.  Merrill,  M.  Schuller,  and  L.  Huang,  15th  STAIF,  (AIP  Press  NY)  1613  (1998). 


5-13 


LIST  OF  SYMBOLS 


Artery  cross-section  area  of  artery  (m2) 

AWall  cross-section  area  of  wall  (m2) 

F  Faraday's  constant  (96,485  C/mole) 

Fij  view  factor  between  surfaces  i  and  j 

hfg  latent  heat  of  condensation  of  sodium  per  unit  mass  (J/kg) 

I  cell  current  (A) 

J  current  density  (A/m2) 

k  thermal  conductivity  (W/mK) 

K  coefficient  includes  conduction  and  radiation  effects  (W/mK) 

ril  sodium  flow  rate  (kg/s) 

P  perimeter  of  cell  wall  (m) 

Pa  pressure  at  anode  side  (Pa) 

Pc  pressure  at  cathode  side  (Pa) 

Pout  electrical  energy  generated  by  the  cell  (W) 

Qair  heat  passes  the  air  as  cooling  fluid  (W) 

Qartery  conducted  heat  to  condenser  through  artery  (W) 

Qcold  heat  loss  at  the  edge  of  cold  plate  (W) 

Qcond  latent  heat  of  condensation  of  sodium  (W) 

Qinsul  heat  losses  to  ambient  through  cell  wall  and  insulation  (W) 
Qloss  total  heat  losses  (W) 

Qrcond  net  radiated  energy  to  condenser  (W) 

Qwall  conducted  heat  to  condenser  through  cell  wall  (W) 

R  perfect  gas  constant  (8314  J/mol.K) 

Rbus  electrical  resistance  of  bus  wire  (Q  m2) 

Rb  ionic  electrical  resistance  (Q,  m2) 

^collector  electrical  resistance  of  current  collector  (£2  m2) 

^contact  contact  resistance  of  the  electrodes  (Q.  m2) 

Rint  total  internal  resistance  (Q.  m2) 

Rop  charge  exchange  polarization  loss  (Q.  m2) 

Rsheet  sheet  resistance  in  the  plane  of  electrode  (Q  m2) 
tfi  thickness  of  BASE  tube  (m) 

T  temperature  (K) 


V  Cell  voltage  (V) 
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cell  voltage  in  open-circuit 
electrode  polarization  over  potential  (V) 
distance  from  condenser  (m) 


emissivity 

cell  conversion  efficiency 
electrical  resistivity  (Q  m) 
electrode  polarization  over  potential  at  anode  (V) 
electrode  polarization  over  potential  at  cathode  (V) 
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Table  1.  Design  Parameters  of  PX-3A  cell 


Cell  diameter  (mm) 

31.75 

Cell  height  (mm) 

101.6 

Evaporator  type 

Deep  Cone 

Evaporator  elevation  (mm) 

5.18 

Evaporator  standoff  thickness 

0.71 

Evaporator  standoff  material 

SS 

Standoff  rings  (mm) 

1.1 

Rings  material 

Ni 

Stud  area  (mm2) 

38 

Stud  material 

SS 

Number  of  BASE  tubes 

5 

Tube  length 

32 

Electrode/tube  (mm2) 

600 

Tube  braze  material 

TiNi 

Current  collector 

60-mesh  Mo 

Feedthrough  braze 

TiCuNi 

Radiation  shield  type 

Circular 

Shield  material 

SS 

Condenser  type 

Creare 

Hot  side 

SS 

Cell  wall 

SS 

Initial  test  date 

7/9/97 

Operation  (hrs) 

12,000 
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Table  2.  Comparison  of  heat  losses  at  the  maximum  power  generated  by  the  AMTEC  cell 


Original 

Case  A 

Case  B 

Case  C 

Case  D 

CaseE 

CaseF 

Pout,  max  (W) 

4.490 

5.971 

5.995 

6.012 

6.109 

6.135 

6.253 

1(A) 

2.447 

3.155 

3.161 

3.166 

3.191 

3.198 

3.228  1 

r\(%) 

12.642 

13.917 

14.278 

14.484 

15.866 

17.067 

20.789 

QinDUt  (W) 

35.514 

42.901 

41.987 

41.511 

38.507 

35.944 

30.078 

Qinsull  (W) 

7.079 

7.492 

7.487 

7.518 

7.099 

7.118 

0.0 

Qarterv  (W) 

3.257 

3.550 

3.207 

3.281 

2.299 

2.291 

2.395 

Qwall  (W) 

5.415 

6.423 

5.674 

5.848 

3.746 

1.115 

1.833 

Qrcond  (W) 

1.300 

1.536 

1.659 

0.860 

1.117 

1.110 

1.246 

Qcond  (W) 

13.495 

17.238 

17.273 

17.298 

17.438 

17.474 

17.643 

Qinsul  (%) 

19.932 

17.465 

17.831 

18.112 

18.436 

19.804 

0.0 

Qarterv  (%) 

9.170 

8.275 

7.639 

7.904 

5.970 

6.374 

7.962 

Qwall  (%) 

15.247 

14.973 

13.951 

14.088 

9.729 

3.103 

6.095 

Qrcon  (%) 

3.662 

3.581 

3.951 

2.072  j 

2.900 

3.088 

4.143 

Qcond  (%) 

38.000 

40.180 

41.139 

41.671 

45.286 

48.615 

58.656 

Table  3.  Comparison  of  maximum  conversion  efficiency  of  the  AMTEC 


Original 

Case  A 

Case  B 

Case  C 

Case  D 

Case  E 

CaseF 

1(A) 

1.98 

2.52 

2.52 

2.44 

2.39 

2.32 

1.93 

^max  (^0 

12.98 

14.38 

14.79 

15.02 

16.64 

18.09 

22.99 

Improvement  in 
efficiency  (%) 

- 

10.79 

13.94 

15.72 

28.20 

39.37 

77.12 

Minimum  BASE 
temperature  (K) 

1062 

1088 

1091 

1093 

1108 

1110 

1126 

Maximum  BASE 
temperature  (K) 

1141 

1153 

1153 

1154 

1156 

1157 

1161 

898 

928 

945 

954 

1010 

1012 

1062 

Maximum  shield 
temperature  (K) 

1003 

1045 

1050 

1054 

1076 

1050 

1103 

Wall  temperature 
at  certain  point  (K) 

1049 

1106 

1110 

mi 

1122 

1138 

1153 
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FIGURE  CAPTIONS 


Figure  1.  A  schematic  diagram  of  vapor-anode  AMTEC  cell. 

Figure  2.  The  electrical  power  generated  by  the  AMTEC  cell. 

Figure  3.  The  heat  supplied  to  the  AMTEC  cell. 

Figure  4.  The  cell  conversion  efficiency. 

Figure  5.  Heat  losses  to  ambient  through  cell  wall  and  insulation. 

Figure  6.  Conducted  heat  to  condenser  through  artery. 

Figure  7.  Conducted  heat  to  condenser  through  cell  wall. 

Figure  8.  Radiation  heat  losses  from  condenser. 

Figure  9.  Latent  heat  of  condensation  of  sodium. 

Figure  10.  Degradation  for  PX-3A  AMTEC  cell  at  different  temperatures. 
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Figure  1.  A  schematic  diagram  of  vapor-anode  AMTEC  cell. 
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Current  (A) 


Figure  2.  The  electrical  power  generated  by  the  AMTEC  cell. 
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Current  (A) 


Figure  3.  The  heat  supplied  to  the  AMTEC  cell. 
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Current  (A) 


Figure  4.  The  cell  conversion  efficiency. 


5-24 


T 


T 


**  *  **  1  *  i  *  >  a  a  a  g  — . .  a  .  a  - 

^  ^  ^  ^  <  <  #-'i<i^.4 . j . j _ j . — 0  0 


.  original  PX3A 

• - •  Case  A 

■ - ■  Case  B 

« - o  Case  C 

a - a  Case  D 

< - «  Case  E 

▼ - ▼  Case  F 


YY — Y-*-Y y — I - y - yJ — y y  '  y — y 1 Y 1 _ ¥-1 _ 1 _ ^ _ i _ Y ' 

1  2  3  4  5  6 

Current  (A) 


Figure  5.  Heat  losses  to  ambient  through  cell  wall  and  insulation. 
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Figure  6.  Conducted  heat  to  cc 
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Figure  7.  Conducted  heat  to  condenser  through  cell  wall. 
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Figure  8.  Radiation  heat  losses  from  condenser. 
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Figure  10.  Degradation  for  PX-3A  AMTEC  cell  at  different  temperatures. 
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SOLUTIONS  OF  THE  TRANSMISSION  LINE  EQUATIONS  USING 
PRODUCT  INTEGRALS  OF  VARIABLE  MATRICES 


Stanly  Steinberg 
Professor  of  Mathematics 
Department  of  Mathematics  and  Statistics 
University  of  New  Mexico 

Abstract 


Product  integrals  and  Lie  algebraic  ideas  are  used  to  find  useful  representations  of  solutions  of  variable 
coefficient  differential  equations  that  have  applications  to  transmission-line  problems. 
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SOLUTIONS  OF  THE  TRANSMISSION  LINE  EQUATIONS  USING 
PRODUCT  INTEGRALS  OF  VARIABLE  MATRICES 


Stanly  Steinberg 


1  Introduction 

The  general  transmission-line  equations  describe  the  currents  and  voltages  on  a  tube  of  transmission  lines: 

=  -mm 

=  -Y{t)V{t)  (1.1) 

where 

t  is  the  position  along  the  tube 
I ( t )  is  the  complex  current  vector  with  N  components 

V  ( t )  is  the  complex  voltage  vector  with  N  components 

Y  ( t )  is  the  N  by  N  complex  admittance  matrix 

Z(t)  is  the  N  by  N  complex  impedance  matrix 


dV(t) 

dt 

dm 

dt 


For  more  details,  see  the  notes  [1]  by  Baum,  Liu,  and  Tesche.  In  the  applications  to  transmission  lines, 
the  vectors  and  matrices  also  depend  on  s,  the  complex  variable  conjugate  to  time  under  the  Laplace 
transform.  Also,  if  we  find  a  complete  set  of  analytic  solutions  then  we  can  easily  solve  problems  with 
inhomogeneous  terms  and  either  initial  or  boundary  conditions,  so  we  do  not  include  these  terms  in  the 
statement  of  the  problem. 

These  transmission-line  equations  (1.1)  can  be  easily  written  as  a  single  second-order  equation  for 
either  the  currents  or  the  voltages: 


dPl  i  i  dl 
__y.y-,__yz/  =  0 


<PV  ,  ,  dV 

^  ~dfi-Z  Z  ~-ZYV  =  ^ 


dt 


(1.2) 


where  Y'  =  dY/dt  and  Z'  =  dZ/dt.  It  is  also  useful  to  write  these  equations  as  a  system  of  2  N  by  2  N 
equations: 


V 

0 

z 

V 

I 

Y 

0 

I 

(1.3) 
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In  fact,  the  matrices  Y  and  Z  cannot  be  arbitrary.  The  reciprocity  condition  requires  that 

YT  =  Y ,  ZT  =  Z.  (1.4) 

To  be  clear,  if  the  complex  conjugate  is  given  by  an  over-bar,  then  the  adjoint  A *  of  a  matrix  A  is  given 
by  A*  =  AT  so  we  do  not  necessarily  have  Y*  =Y  or  Z*  =  Z .  So  in  the  general  case  we  do  not  assume 
the  matrices  are  Hermitian  (self-adjoint).  The  condition  that  there  are  no  infinite  physical  parameters 
requires  that  the  Y  and  Z  matrices  be  invertible,  that  is, 

Y~l  and  Z~l  exist.  (1.5) 

As  noted  above,  the  matrices  Y  and  Z  are  functions  of  a  parameter  s.  For  apropriate  ranges  of  s,  in 
particular  for  s  purely  imaginary,  these  matrices  have  positive  real  part,  so  we  just  look  at  that  case  and 
assume  that  for  all  vectors  V  find  I 

ft  (VtY V)  >  0  and  ft  (/T ZI)  >  0 .  (1.6) 

To  summarize,  we  will  always  assume  that  the  matrices  in  the  (1.1)  are  symmetric  and  have  strictly 
positive  real  part. 

There  are  other  special  cases  of  interest.  When  the  materials  in  the  cables  are  uniform,  then  the 
matrices  are  constant  in  t, 

Y{t)=Y  and  Z(t)  =  Z .  (1.7) 

In  the  lossless  cables  case,  the  vectors  and  matrices  are  real  for  real  s  and  the  solutions  for  all  s  can  be 
found  from  s  =  i  by  rescaling,  so  we  can  take 

Y{t)=iL(t)  and  Z(t)  =  iC{t).  (1.8) 

where  L  is  the  inductance  matrix  and  C  is  the  capacitance  matrix  and  both  L  and  C  are  symmetric  and 
thus  Hermitian  (self-adjoint)  and  are  positive.  The  most  specific  case  of  interest  is  when  there  is  uniform 
permeability  ji  and  uniform  permittivity  e  in  which  case 

L  =  fiF  and  C  =  eF~l  (1.9) 

where  F  only  depends  on  the  geometry  of  the  cables  and  is  real  symmetric,  so  Hermitian  and  positive 
definite. 
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1.1  Uniform  Cables 


In  the  case  of  uniform  cables,  we  can  write  the  solution  of  the  transmission-line  equations  (1.1)  using  either 
hyperbolic  or  trigonometric  functions.  The  solution  of  (1.1)  is  given  by  an  exponential  of  the  coefficient 
matrix  found  in  (1.3): 


V(t)  I  Y  0  I  V(0) 

=  e  L  J  .  (1.10) 

.  *(*)  J  L  /(°)  . 

Now  if  we  write  the  exponential  as  a  power  series  and  collect  every  other  term  into  two  groups,  we  obtain 
0  Z 

-t  r  _  _  _ 

Y  0  cosh  (tVZY)  -y-1v/ZFsinh(tVzF) 

e  L  J  =  _  .  (1.11) 

—Z  1  y/ZY  sinh(t  VZY)  cosh  (t\/ZY) 

If  we  write  Y  =  iL  and  Z  =  i  C,  where  L  is  the  inductance  matrix  and  C  is  the  capacitance  matrix,  and 
both  L  and  C  are  symmetric  and  thus  Hermitian  (self-adjoint),  then 


[  L  0  J  _  cos(ty/CL)  -iL-'s/CLsinitVCL) 

—i  C-1  \JCL  sin(t  V CL)  cos  (t  VCL) 

So,  in  this  case,  the  solution  is  given  by  simple  waves  traveling  down  the  line. 


2  Symmetrization  of  the  Transmission  Line  Equations 

In  some  cases,  the  Equations  (1.3)  can  be  symmetrized  by  a  change  of  dependent  variables.  To  simplify, 
assume  the  material  properties  aire  constant  and  then  consider 


d  AV 

dt  d  t 


0  AZB~ 


BYA~ 


where  the  matrices  A  and  B  axe  to  be  chosen  so  that  the  matrix  in  the  previous  equation  is  symmetric, 
that  is,  so  that 

AZB~1=BYA~1,  (2.14) 

or 

Y  =  B~1AZB~1A.  (2.15) 

Setting  M  =  B_1A,  given  Y  and  Z  we  must  find  M  so  that 


Y  =  MZM. 
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This  is  a  system  of  N 2  quadratic  equations  in  N 2  unknowns,  so  it  is  not  clear  if  there  even  exist  complex 
solutions  or  not.  However,  if  Y  and  Z  commute  then  there  is  a  solution: 

[Y,Z)  =  Y  Z -ZY  =  0^  M  =  \/Y\/Z~1  =VZ~1  VY .  (2.17) 

So  in  the  case  of  commuting  matrices,  the  simple  choice  of 

A  =  VY,  B  =  VZ  (2.18) 


will  make  the  system  symmetric. 


2.1  Variable  Coefficients  and  Symmetry 

We  now  assume  that  for  each  z,  the  matrices  in  the  transmission  line  equations  (1.1)  commute: 

[m  Z(t))  =  Y(t)  Z(t)  -  Z(t)  Y(t)  =  0 .  (2.19) 

If  we  introduce  the  new  variables 


V(t)  =  A(t)V(t),  I(t)  =  B(t)I(t), 


(2.20) 


then  a  simple  computation  gives 

dV{t) 

dt 

d/(t) 

dt 


^  A~l  (t)  V  -  A(t)  Z(t)  B~x  (t)  /«) 
^  B -1  (i)  /  -  B(t)  Y(t)  A"1  (t)  V(t) 


These  equations  can  be  written  as  a  system 


d 

V 

i 

A{t)  Z(t)  B~l(t) 

V 

dt 

i 

.  B(t)Y(t)A~'(t) 

B'(t)B~'(t) 

I 

The  choice  of 


A(t)  =  ^/m,  B{t)  =  ^/m 


(2.21) 


(2.22) 

(2.23) 


makes  the  off-diagonal  part  of  the  system  symmetric,  and  the  diagonal  parts  are  logarithmic  derivatives. 


2.2  Problems  With  Non-Commuting  Matrices 

Here  we  indicate  why  the  dealing  with  matrix  functions  of  t  that  do  not  commute  is  difficult.  So,  if  a 
matrix  function  M(t)  has  the  property  that  the  matrices  for  different  values  of  t  commute,  that  is, 

[M(ty),M(t2)}  =  0  (2.24) 
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for  all  ti  and  t2 ,  and  f(x)  is  any  differentiable  function,  then 

jtf{M{t))  =  £(MW)  itr{t)  ’  (2-25) 

that  is,  the  chain  rule  works.  Without  the  commuting  condition,  the  chain  rule  does  not  work  in  general, 
as  elementary  examples  will  show. 

There  is  one  important  exception  for  the  function  1/x.  If  we  differentiate  M{t)  M~x{z)  =  Identity, 
the  we  get 

jM-1®  =  — M~\t )  ( jtM{t ))  (2.26) 

One  important  problem  is  that  the  derivative  of  a  square  root  of  a  matrix  cannot  be  found  in  closed 
form,  but  we  can  find  a  simple  equation  for  the  square  root.  To  see  this,  we  follow  the  calculus  derivation 
of  the  derivative.  We  have  that 

M(t)  =  y/Y(f)^  M2{t)  =  Y(t).  (2.27) 


and  then  differentiation  gives 


dM(t) 

dt 


M(t)  =  Y'(t). 


(2.28) 


This  doesn’t  lead  to  a  simple  formula  for  the  derivative  of  the  square  root.  It  does  give  a  N  x  N  system 
of  nonlinear  equations  to  find  the  derivative  of  the  the  square  root,  but  again  it  is  difficult  to  know  when 
this  system  is  solvable. 

We  would  also  like  to  have  a  simple  form  for  the  “logarithmic”  derivatives 


dA(t)  dB(t)  ^ 

dt  ’  dt 


(2.29) 


that  appear  in  the  transformed  equations.  However,  these  are  not  the  derivatives  of  ln(j4(£))  and  ln(B(t)) 
unless  the  commuting  condition  (2.24)  holds.  If  we  set 


M(t)=ln(A(t)), 


(2.30) 


then 

eM(t)  _  _  (2.31) 

This  last  expression  can  be  differentiated  using  Formula  6.196  in  [4],  but  this  produces  a  formula  difficult 
to  solve  for  dM /dt. 
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2.3  Product  Integrals  and  Logarithmic  Derivatives 

The  theory  of  product  integrals  and  logarithmic  derivatives  was  invented  to  deal  with  just  the  situation 
described  in  the  previous  section.  Using  the  definition  of  the  product  integral  (see  Section  5  of  Baum  and 
Steinberg  [4]  for  definitions  and  basic  properties)  applied  to  the  exponential  of  Af(f),  we  define 

t 

G(t)=ne<,*MW>  (2-32) 

o 

and  then  we  get  (see  [4],  Equation  (5.155)) 

^(t)  =  M(t)G(t).  (2.33) 

Also  (see  [4],  Equation  (5.158))  the  logarithmic  derivative  is  defined  as 


DtG(t)  =  G'(t)G-'(t), 


(2.34) 


so  that 

DtG(t)  =  G'(t)G-1(t)  =  M(t). 
We  also  have  that  (see  [4],  Equation  (5.160)) 

t 

G(t)  =  JJ  eD*G(g)  ds . 


(2.35) 


(2.36) 


One  of  the  more  interesting  formulas  for  product  integrals  is  the  analog  of  the  formula  exp  (a  +  b)  = 
exp  (a)  exp  (5)  ((see  [4],  Equation  (5.163)):  if 

t 

(2.37) 


G(t)  =  II eMt)  d"  ’ 


then 


jQ  eM(*)+B(*))  *  =  J"JeG  G(«)  ds 


(2.38) 


We  will  use  this  formula  in  the  next  section. 


2.4  Product  Integral  Reduction  of  the  Transmission  Line  Equations 


Here  we  show  that  the  sum  rule  (2.38)  can  be  used  to  check  that  we  have  correctly  transformed  the 
transmission  line  equations  or,  conversely,  that  the  sum  rule  is  correct.  We  rewrite  Equation  (2.22)  as 


V 

( 

’  A' A~l 

0 

0 

AZB-1 

\ 

V 

= 

+ 

i 

V 

0 

B'B”1 

BY  A-1 

0 

) 

} 

(2.39) 
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so  and  then  apply  the  sum  rule  (2.38)  with 


so  that 


A' A-1  0 

0  B'B-1 

0  AZB- 
BYA-1  0 


G  -4 


G~lBG 


Y  0 


and  then 


A,{s)A-1(s)  .A(s)  Z(s)B-1(s) 

B'(s)B-1(s) 


A(t)  0  ‘ 

11  e 

0  B(t)  o 


0  Z(s) 
Y(s )  0 


,  (2.42) 


which  agrees  with  the  definition  (2.20)  which  we  can  write  as 

v(t)  i  =  r  A(t)  o  l  r  v(t ) 

_  i(t)  I  0  B(t)  I(t) 


3  Diagonzdizing  Symmetric  Transmission  Line  Equations 

In  this  section  we  assume  the  commuting  property  (2.19) 

[m  z(t)]  -  Y(t)  Z(t)  -  Z(t)  Y(t)  =  0 .  (3.44) 

and  that  one  of  the  matrices  Y  (t)  or  Z(t)  can  be  diagonalized.  Then,  because  of  the  commuting  property, 
both  of  the  matrices  can  be  diagonalized  and  have  a  complete  common  set  of  eigen- vectors  for  all  t. 

Define  the  S  matrix  to  be  the  inverse  of  the  matrix  with  columns  the  eigenvectors: 

S(t)  =  [&(*),•• -Mt)]"1,  (3.45) 


and  then  we  have 


A r(t)  =  S(t)  Y(t)  S-'lt) ,  A z(t)  =  5(f)  Z(t)  S'^t) , 
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where  Ay  and  A z  are  diagonal.  If  we  now  introduce  the  variables 


i(t)  =  smt),  m  =  s(t)v(t), 


(3.47) 


then  the  transmission  line  equations  can  be  written  as  in  (2.22): 

£  [  V  1  _  [  S'iQS-'it)  5(t)Z(t)5-J(t)  j  [  V  " 

dt[  I  \  [  5(t)y(t)5-x(t)  5'(t)5-x(t)  J  [  _ 

The  eigen-decomposition  (3.46)  then  gives 


5'(t)5-x(t)  A z(t)  1  [  V  ' 

A Y{t)  S'(t)  S-'it)  J  [ 


If  the  diagonal  logarithmic  derivative  terms  are  zero,  then  the  system  is  uncoupled: 


(3.48) 


(3.49) 


£  r  v  j  _  r  a z(t)  j  r  v 

dt[~I  \  [  AY(t)  0  J  [ 


(3.50) 


These  terms  will  be  zero  if  the  eigenvectors  are  constant.  If  they  are  not  zero,  then  the  system  doesn’t 
uncouple  and  the  standard  solution  techniques  cannot  be  applied. 


3.1  Circulant  Matrices 

Circulant  matrices  give  an  ideal  situation  in  which  to  apply  the  diagonalization  techniques  because  they 
have  constant  eigenvectors.  Circulant  matrices,  which  are  constant  along  “diagonals”,  are  discussed  in 
detail  in  Nitsch,  Baum  and  Sturm  [3]  and  Davis  [2].  Here  we  adopt  a  definition  in  terms  a  basis  which  is 
convenient  for  Lie  algebraic  Computations.  In  particular,  we  will  define  circulant  matrices  in  terms  of  the 
shift  operation. 

As  shown  in  Figure  3.1,  put  N  evenly  distributed  points  labeled  by  1, 2,  •  •  •  N  —  1,  N  on  the  unit  circle 
and  then  assigning  the  values  f  =  {/i,  /2,  •  •  ■  In}  to  these  points.  The  shift  S  =  S(N)  in  the  clockwise 
direction  moves  the  value  at  the  point  1  to  the  point  N  —  1,  the  value  at  the  point  2  to  the  point  1  and 
so  forth.  So  if  g  =  Sf,  then 


5x  =  (Sf) i  =  h  , 
92  =  (Sf)2  =  h , 

9n  =  (S/)n  =  fi  ■ 


(3.51) 

(3.52) 
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In  fact,  we  will  define  a  family  of  shift  operators  Sh(fc)  =  Sh(fc,  N): 

{1  if  mod(i  -  j  +  k,  N)  =  0 

V  '  (3.56) 

0  if  mod(i  —  j  +  k,  N)  ^  0 

The  geometric  interpretation  of  Sh(fc)  is  that  the  values  at  the  points  on  the  unit  circle  are  shifted 
clockwise  k  points  if  k.geqO  and  counter  clockwise  —k  points  if  k  <  0.  The  matrices  Sh(k,N)  are  distinct 
for  0  <  k  <  N  and  Sh(fc  +  mN,N)  =  Sh(k,N)  for  all  m.  Also  Sh(0,JV)  is  the  identity  matrix  and 
Sh(l,IV)  =  Sh(l)  =  S(N)  =  S. 

For  N  =  5  the  distinct  shift  operators  are: 


Sh(0) 


Sh(2)  = 


1 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

1 

0 

0 

,  Sh(l)  = 

0 

0 

0 

1 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

1 

1 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

1 

,  Sh(3)  = 

1 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

1 

0 

0 

0  0  0  0  1 

1  0  0  0  0 

Sh(4)=  0  10  0  0 

0  0  10  0 
0  0  0  1  0 

From  the  geometric  interpretation  of  the  shift  operator,  it  is  clear  that 

Sh{kuN)Sh{k2,N)  =  Sh(kl+k2,N). 


(3.57) 


(3.58) 


(3.59) 


(3.60) 


An  algebraic  proof  of  this  requires  only  a  little  work  with  the  mod  function.  This  immediately  implies 
that 

Sh(fc,  N)  =  Sk(N) ,  (3.61) 

for  all  integers  k. 
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(3.62) 


A  circulant  matrix  is  a  matrix  that  can  be  written  as  a  linear  combination  of  shift  matrices: 

JV-l 

Ci(co,ci,  •  •  •  cjv-i)  =  ^2  c*  Sk  • 
k= o 

Such  matrices  are  constant  on  along  the  “diagonals”.  It  is  show  in  [2]  that  other  common  definitions  of 
circulant  are  equivalent  to  this  definition.  In  any  case,  all  five  dimensional  circulant  matrices  are  given  by 

Ci(co,ci,C2,C3,c4)  =  co  Sh(0)  +  ci  Sh(l)  +  c2  Sh(2)  +  c3  Sh(3)  +  c4  Sh(4)  (3.63) 

Co  Cl  C2  C3  C4 

C4  Co  Ci  C2  C3 

C3  C4  Co  Ci  C2  •  (3.64) 

C2  C3  C4  Co  Cl 

Ci  C2  C3  C4  Co 

3.2  Diagonalizing  Circulant  Matrices 

Because  circulant  matrices  can  be  written  as  a  polynomial  in  the  shift  S  (3.62),  they  all  commute  under 
multiplication  and  consequently  can  be  diagonalized  simultaneously.  If  fact,  if  we  find  a  basis  that  diago¬ 
nalizes  S,  then  the  same  basis  diagonalizes  all  powers  of  S,  and  consequently  all  circulant  matrices.  It  is 
well  known  that  the  discrete  Fourier  transform  diagonalizes  all  translation  invariant  operators,  that  is,  all 
operators  the  commute  with  the  translation  S,  so  lets  look  at  the  discrete  Fourier  transform 
Let  u  be  a  primitive  JV-th  root  of  unity.  Typically  we  take 
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There  are  many  choices  possible  for  the  discrete  Fourier  transform  matrix.  Our  choice  gives 

F*S F  =  diag(l,  w~l ,  w~2,w~3,  •  •  • ,  wx^N) ,  (3.70) 

that  is,  the  fc-th  column  of  F  is  an  eigenvector  of  S  with  eigenvalue  ul~k .  Additionally,  F  is  a  unitary 
matrix. 

When  we  solve  the  transmission-line  equations  using  Lie  techniques,  a  critical  ingredient  is  the  expo¬ 
nential  of  each  basis  element.  For  circulant  matrices,  a  basis  is  S*,  0  <  k  <  N.  We  have  diagonalized  5, 
so 

F*Sk  F  =  diag(l,  w~k,  w~2  k,  w~3  fe,  •  •  • ,  k) ,  (3.71) 

and  then  we  see  that 


2  It  *.,,-3  1. 


e* =  F  diag(l,  e* »  ' ,  e* ' u,~* " ,  e* ’ e* 1 u,<1~'v) 1  )F* . 


(3.72) 


To  apply  this  to  the  transmission-line  equations  (1.3), 


dt 


we  need  to  write  Z  and  Y  in  terms  of  S: 

TV— 1 


V 

0 

z 

V 

I 

Y 

0 

I 

TV-1 


Z(t)=J2h(t)Sk,  Y(t)=£ok(t)S*. 


k= 0  k= 0 

In  this  case  the  system  of  equations  can  be  rewritten  as 


d_ 

dt 


V 

I 


TV-1 

=  X  -«*(<)  W) 

i,j=0 


o 

1 _ 

V 

|_  S  *  0 

I 

so  we  only  need  to  know  how  to  exponentiate  the  matrices 

0  S> 

,  0  <  i,j  <  N  -  1 . 

S *  0 

However, 


F* 

0 

0  Sj 

F  0 

0 

(F*SF)i 

0 

F* 

o 

_ i 

0  F 

(F*SF)* 

0 

(3.73) 


(3.74) 


(3.75) 


(3.76) 


(3.77) 


which  gives  a  “skew  diagonalization”  of  the  matrix  involved  in  the  exponential. 

For  the  “skew-diagonal”  types  of  matrices  involved  in  the  transmission-line  equations,  the  xxx 
Lemma: 


7-14 


4  Conclusions 
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APPROXIMATION  OF  THE  SPACE-CHARGE  LIMITED  CURRENT  OF  RELATIVISTIC 
ELECTRON  BEAMS  IN  COAXIAL  DRIFT  TUBES 


Kenny  F.  Stephens  II 
Doctoral  Candidate 
Department  of  Physics 
University  of  North  Texas 


Abstract 

A  two-dimensional  theory  for  predicting  the  space-charge  limited  current  of  relativistic  electron 
beams  in  coaxial  drift  tubes  is  presented.  Applicable  to  annular  beams  or  arbitrary  radius  and 
thickness,  the  theory  considers  not  only  the  effect  of  the  inner  conductor  but  also  takes  into 
account  finite-length  effects.  Using  Green’s  second  identity,  Poisson’s  equation  for  the  space- 
charge  limited  problem  is  reduced  to  solving  a  Sturm-Liouville  eigenvalue  problem.  For  some 
limiting  cases,  closed-form  analytical  solutions  are  available. 
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APPROXIMATION  OF  THE  SPACE-CHARGE  LIMITED  CURRENT  OF  RELATIVISTIC 
ELECTRON  BEAMS  IN  COAXIAL  DRIFT  TUBES 

Kenny  F.  Stephens  II 


I  Introduction 


One  limitation  to  the  production  of  high-power  microwaves  (HPM)  is  the  ability  of  the  HPM 
device  to  carry  an  electron  beam  of  sufficient  current.  It  is  well  known  that  the  space-charge  of  the 
beam  places  an  upper-limit  on  the  amount  of  current  than  can  be  carried.  As  the  beam  current 
approaches  this  limit,  the  amount  of  current  actually  propagated  through  the  device  no  longer 
increases.  Eventually,  only  a  fraction  of  the  input  current  travels  through  the  device  regardless  of 
the  amount  of  input  current.  Because  this  phenomena  depends  on  the  space-charge  of  the  beam, 
this  limiting  value  is  called  the  space-charge  limited  (SCL)  current.  Thus,  the  ability  to  predict 
the  SCL  current  is  important  in  a  variety  of  intense  electron  beam  applications,  including  HPM 
devices. 

Precisely  solving  the  SCL  problem  requires  the  solution  of  a  highly  nonlinear  Poisson  equation. 
To  date,  only  one-dimensional  solutions  are  available  in  closed,  analytic  form.  Because  of  these 
limited  solutions,  the  researcher  cannot  consider  both  radial  and  axial  effects  simultaneously 
without  resorting  to  numerical  solution  or  computer  simulations.  For  this  reason,  attempts  to 
estimate  the  SCL  current  are  important.  Perhaps  the  simplest  approach  is  the  “infinite-length” 
approximation.  Under  this  approximation,  the  beam  is  considered  to  travel  in  a  smooth-bore  bore 
under  the  influence  of  an  infinitely  strong  axial  magnetic  field  that  prevents  radial  and  azimuthal 
motions.  Furthermore,  all  physical  quantites  are  considered  independent  of  the  axial  coordinate  so 
that  only  radial  variations  in  the  electric  potential  need  be  considered.  This  approach  is  adequate 
for  situations  where  the  beam  length  is  several  times  the  tube  radius. 

In  Ref.  [1],  a  two-dimensional  theory  is  presented.  Assuming  the  beam  travels  in  a  strong  axial 
magnetic  field  so  that  radial  and  azimuthal  motions  need  not  be  considered,  the  theory  presents 
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a  method  to  estimate  an  upper-bound  on  the  SCL  current.  Applicable  only  to  beams  that 
completely  fill  a  hollow  drift  tube,  the  theory  was  later  extended  to  annular  beams  of  arbitrary 
radius  and  thickness  in  Ref.  [2],  Since  the  use  of  coaxial  drift  tubes  has  recently  been  shown 
to  enhance  the  SCL  current  (see  Ref.  [3]),  the  present  work  extends  the  theory  to  consider  the 
finite-length  effects  of  relativistic  electron  beams  in  coaxial  drift  tubes. 


II  Theory  Review 


The  problem  geometry  consists  of  two  coaxial,  cylindrical  conductors  of  length  L.  The  inner 
conductor  has  radius  rq  and  the  outer  conductor  radius  r<i-  The  annular  electron  beam,  of  inner 
radius  ?q,  outer  radius  r0  and  current  Jo,  is  emitted  from  a  foil  at  2  =  0  with  uniform  density 
and  passes  to  a  completely  absorbing  collector  plate  at  z  =  L.  A  strong,  axial  magnetic  field  is 
externally  applied  to  prevent  any  radial  or  azimuthal  motion.  Letting  the  beam  be  injected  with  a 
kinetic  energy  [70  —  l]mc2,  the  kinetic  energy  at  any  location  within  the  drift  tube,  [7(7",  z)  —  l]mc2, 
can  be  written  in  terms  of  the  injected  kinetic  energy  and  the  electrostatic  potential,  z), 
through  conservation  of  energy: 


7  (r,z)  =  70  +  — ~<t>{ryz) 


me * 


(i) 


Due  to  the  external  magnetic  field,  the  relativistic  mass  factor  is  simply  related  to  the  axial 
velocity  according  to 


7  (r,  z)  = 


(2) 


where  /3c  «  vz  is  the  axial  velocity.  Since  only  axial  motion  is  considered,  the  current  density  is 
simply  J  =  Jz  —  —  e7io/?c,  where  no  is  the  beam  density.  Through  conservation  of  charge,  the 
charge  density  must  be  constant.  Using  Eqs.  (1)  and  (2),  Poissong  equations,  V2^(r,  z)  =  47reno, 
can  be  written 

,  n  <  r  <  r0, 0  <  z  <  L, 

=  <  rv,z)  (3) 


V  *(r,z)  =  -^;r 


\  d  d<b  (r,  z)  d2§  (r,  z ) 


dr 


dz 2 


P  (r,z)’ 

0,  otherwise 
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subject  to  the  boundary  conditions 

$  (r,  0)  =  $  (r,  L)  =  0,  r\  <r  <r^ 

(4) 

$  (ri,  z)  =  $  (r2,  z)  =  0,  0  <  2  <  L. 

Here,  $  =  etp/mc?  is  the  normalized  potential  energy  and  K  —  4io/(r^  —  rf)  with  io  =  Io/ 1  a  the 
injected  current  normalized  to  the  Alfven  current,  I  a  —  mc3/e.  The  velocity  factor  in  Poisson’s 
equation  is  written  in  terms  of  the  potential  according  to 


f3(r,z)  =  Jl  -7  2{r,z)  = 


1  ~  [7o  +  £  (r,  z)Y 

7o  +  $  (r,  2) 


Equation  (3)  is  the  nonlinear  system  whose  solution  is  desired.  Although  an  exact  solution  is 
unavailable,  the  next  section  reviews  the  method  developed  in  Ref.  [1]  which  will  yield  an  upper- 
bound  on  the  SCL  current. 


Ill  Upper  Bound  for  the  SCL  Current 


To  make  progress  toward  an  approximate  solution  of  Eq.  (3),  Green’s  second  identity  is  applied 


to  the  beam  volume: 


[  Uv2$  -  $V2V>1  dv  =  /  lip—  -  ds . 

Jv  L  J  Js  L  on  on_ 


The  underbraces  indicate  that  r  represents  the  volume  integral  and  E  the  surface  integral.  In 
Eq.  (6),  is  the  electrostatic  potential  defined  in  Eqs.  (3)  and  (4)  while  ^  is  a  solution  to  the 
eigenvalue  problem  defined  by 


V  ^(r,z)  =  -—r 


1  d  cty>(r,  z)  d2/ip(r,z) 


r  dr  dr 


AV>  (r,  z) ,  ri  <  r  <  r0j  0  <  z  <  L,  (7) 


ip  (r,  0)  =  ip  (r,  L) 


dty  (r,  z)  i  f  \ 

am  — ^ -  -  a2ip  z) 

or  r—Ti 

hro^A  _w(ro,2) 

(7r  r=r0 


=  0,  n  <r  <  r0. 
=  0,  0  <  z  <  L, 


=  0,  0  <Z<L. 
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The  constants  a  and  b  are  determined  in  Sect.  IV  such  that  1)  the  surface  integral  E  in  Eq.  (6) 
is  zero,  2)  Eq.  (7)  has  an  eigenfunction  that  is  nonzero  within  the  problem  domain,  and  3)  the 
corresponding  eigenvalue  is  negative.  Assuming  these  conditions  are  satisfied,  the  volume  integral 
in  Eq.  (6)  must  also  be  zero  and  can  be  written  as 


27T 


j  r 

Jo  J n 


K 


A0<J>  (r,  z) 


ipo  (r,  z)  rdrdz  =  0, 


(9) 


[0  (r,  z ) 

where  Ao  is  the  fundamental  eigenvalue  corresponding  to  the  eigenfunction  rpo(r,z)  that  is  zero 
only  on  the  volume  surface.  Since  ipo  is  does  not  change  sign  within  the  integration  volume,  the 
bracketed  term  in  Eq.  (9)  must  change  sign.  Let  the  point  at  which  this  occurs  be  ( r',z Thus, 


K  =  \^\0{r\z')§{r\z'). 


(10) 


If  the  derivative  of  K  with  respect  to  7 (r',z')  is  set  equal  to  zero  and  solved  for  7 (r\zf),  it  is 
found  that  /3(r\  z')  is  maximal  when  7 (r\zf)  =  (7q//3  —  l)3/2.  Since  1  <  7 (r\zf)  <  70, 

this  implies  that  K  has  the  upper  bound 


^  <|A0|  (7o2/3-l)3/2. 


(11) 


Therefore,  Eqs.  (10)  and  (11)  define  the  upper  bound  for  the  SCL  current  as 

-2  *  -.3/2 


r2 

io  <iub  =  — 


-|Ao|  (ll,Z  -  V 


(12) 


Thus,  the  problem  of  determining  an  expression  for  the  SCL  current  for  a  specified  beam  and 
drift  tube  geometries  is  reduced  from  solving  a  nonlinear  partial  differential  equation  to  solving 
a  two-dimensional  Sturm-Liouville  eigenvalue  problem.  The  general  solution  of  this  problem  is 
developed  in  the  following  section. 


IV  General  Solution  of  the  Eigenvalue  Problem 


An  upper  bound  for  the  SCL  current  was  obtained  in  Sect.  Ill  based  on  the  assumptiong  that 
the  surface  integral  in  Eq.  (6)  is  zero.  This  section  solves  the  eigenvalue  problem  defined  by  Eqs. 
(7)  and  (9)  such  that  the  surface  integral  is  indeed  zero.  This  is  achieved  by  determining  the 


8-6 


coefficients  ai,  0,2,  b\  and  62-  By  separation  of  variables,  the  desired  eigenfunction  can  be  written 


ip(r,z)  =  p(r)  sin  (kiz)  ,  (13) 

where  Kn  —  mr/L.  This  form  automatically  satisfies  the  axial  boundary  conditions.  Substituting 
Eq.  (13)  into  E  yields 


E  =  27t£  [L  p(t) 
Jo 


dr  |r=f 


sin  {k\z)  dz\ 


The  solution  to  Eq.  (3)  in  the  charge-free  regions  as 


where 


<t>  (r,  z)  =  ]T  An9n  ( r ,  ri)  sin  (Knz) ,  rx  <  r  <  n, 

n— 1 
oo 

<P  (r>  Z)  =  Y1  Bn9n  {r,  r2)  sin  (Knz) ,  ra  <  r  <  r2, 


9n  (r,  r')  =  I0  ( Knr )  -  ^~~rK0  (Knr) . 

\^nr  ) 


Io  and  Ko  are  the  modified  Bessel  functions  of  the  first  and  second  kind,  respectively.  The  form 
of  the  function  gn(r)rf)  is  chosen  to  satisfy  the  radial  boundary  conditions  and  to  reduce  to  the 
proper  form  for  hollow  tubes,  r\  =  0. 


If  Eq.  (15)  is  substitued  into  Eq.  (14),  it  is  found  that  the  surface  integral  does  vanish  is  ip 
non-trivially  solves  the  eigenvalue  problem 

d  LdP(r)l  ,  .  A 


—  r  —  +  Arp  (r)  =  0,  n  <  r  <  r0, 


amp'  {n)  -  a2p  ( n )  =  0, 


hr  op'  (r0)  -  b2p(r0)  =  0. 

where  the  coefficients  in  the  boundary  conditions  are  given  by 


ai  =  Si  (n,ri) 

dg\  (r,  rx)  I 


02  =  T{ 


h  =  gi  (r0,  r2) 

te  =  ro*!jpi2) 


Noting  that  the  differential  equation  for  p(r)  is  Bessel’s  equation,  the  solution  to  Eq.  (17)  can  be 
written 


P  (r)  =  Jo  (kr)  +  CN0  (hr) , 
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where  Ao  =  k 2  and  Jo  (No)  is  the  zeroth  order  Bessel  functions  of  the  first  (second)  kind.  To 
determine  the  coefficient  C  and  eigenvalue  k,  the  following  system  of  equations  must  be  solved: 


a2  [Jo  (kri)  +  CN0  (kri)}  +  aj/crj  pi  (kri)  +  CN\  (kri)]  =  0, 
b2  [Jo  (kr0)  +  CN0  ( kr0 )]  +  bikra  [Ji  ( kr0 )  +  CN\  ( kra )]  =  0. 


(21) 


In  some  cases  it  is  not  necessary  to  determine  the  value  of  C  so  that  only  the  equation 

a2Jo(kri)  +  aikriJi(kri)  _  b2  Jo  (kr0)  +  bikr0J\  (kr0) 
a2N0  (kn)  +  aik^Ni  (kri)  b2N0  (kr0)  +  bikr0Ni  (kra)  ’ 

need  be  solved  for  k. 


(22) 


Having  now  solved  the  eigenvalue  problem  defined  in  Sect.  Ill,  the  upper  bound  for  the  SCL 
current  can  be  written  directly  in  terms  of  the  eigenvalue  k  as 

iuB  =  —  4  ~  [fc2  +  k2]  (to/3  -  i)  /  •  (23) 

To  estimate  the  SCL  current  for  an  annular  electron  beam  travelling  through  a  coaxial  drift  tubej 
the  value  of  k  must  be  determined  from  Eq,  (22)  and  inserted  into  Eq.  (23).  The  next  section 
presents  some  closed-form  analytical  expressions  for  the  upper  bound  under  certain  limiting  cases. 


V  Coaxial  drift  tubes 


Having  developed  an  expression  for  the  SCL  current  upper  bound  in  terms  of  the  eigenvalue 
problem  discussed  in  Sect.  IV,  the  determinant al  equation  for  the  eigenvalue  is  now  solved  for 
cases  which  yield  closed  analytic  expressions.  This  section  determins  the  SCL  current  for  annular 
beams  in  coaxial  tubes.  For  cases  of  hollow  tubes,  the  reader  is  referred  to  Ref.  [2]. 


V.l  Thin  annular  beam 

The  first  case  is  a  thin  annular  beam,  i.e.,  a  beam  with  zero  thickness.  Letting  the  beam  thickness 
be  5  —  r0  —  r*,  the  inner  beam  radius  is  replaced  with  r;  — >  (1  —  8)r0  in  the  boundary  condition 
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coefficients  a\  and  02-  These  coefficients  are  then  replaced  with  their  first-order  Taylor  series,  i.e., 

a-i  ~  9x  {r0,  n)  ~  KiroSg'i  (r0,  rx) , 


where 


02  ~  «i r0  [g'l  (r0,r  1)  -  n\r08gi  (r0,ri)]  , 


Mr»,n  )=^ 


Next,  the  circular  Bessel  functions  in  Eq.  (22)  are  replaced  by  the  leading  term  of  their  respective 
asymptotic  expansions,  yielding 


cot (kS)  —  — 


a2&2  +  gjji k2r2  (1  -  5) 
a\b2k  (1  -  5)r0  —  a2bikr0  ’ 


Considering  only  the  right-hand  side,  replaced  a\  and  a2  Using  Eq.  (24,  the  right-hand  side  can 


be  re-written  as 


cot ( k5 )  = 


_ b\r0  [Jq  («i r0)  -  K0  (Kir0)] _ 

binlr0  [Ii  («i rQ)  +  Ki  («irc)]  b2  [J0  («i rQ)  +  K0  (*ir0)]  ’ 


where  the  limit  5  — >  0  was  taken  and  only  the  leading  term  was  retained.  Expanding  b\  and  b2: 
the  modified  Bessel  function  Wronskian 


I1(z)Ko(z)  +  I0(z)Kl  (z)  = 


allows  the  determinantal  equation  to  be  written 

cot  (M)  =  («iri)l  (W) .  (29) 

L  o  gi{ri,r2)  J 

To  simply  further,  note  that  the  right-hand  side  is  a  linear  function  of  (kS)  whose  slope  increases 
without  bound  as  S  approaches  zero.  Thus,  equality  will  hold  only  if  kS  <C  1.  Using  the  small 
argument  expansion  of  the  cotangent,  cot (z)  ~  l/z,  gives 


k2Sr0 


_ g\  (nyrg) _ 

9i  (' r0 ,  ri)  g\  (r0,  r2)  K0  (^m) ' 


Having  obtained  a  closed  expression  for  A:2,  the  upper  bound  SCL  current  can  now  be  determined. 
If  the  beam  area  is  re-written  as 


7T  (r20  -  r?)  =7 r  (r0  -  r4)  (rQ  +  r4)  =  25r0, 
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the  geometric  factor  is 


9i  (ri,r2) 


9\  {r0,  n)gi  (r0,  r2)  K0  («i n) 


and  the  SCL  current  estimate  is 


WB  = 


9i  (ri,r2) 


9i  (■ rb,ri)gi  (r6)r2)/<'o(«m) 


(7o2/3-l)3/2. 


(32) 


(33) 


V.2  Thin  annular  beam,  long  drift  tube 


The  result  in  Eq.  (33)  is  now  shown  to  agree  with  the  infinite-length  approximation.  To  this 
end,  the  beam  is  further  assumed  to  be  long.  This  implies  that  the  beam  length  is  larger  than 
any  radial  dimension,  i.e.,  r 2  L.  With  the  small  argument  expansion  of  the  modified  Bessel 
functions, 

Io(z)  ~  1, 

K0  (z) - In  (z) , 

the  geoemtric  factor,  Eq.  (32),  becomes 

*"fe) 


(34) 


21»fe)tofe)' 

Therefore,  the  SCL  current  estimate  for  a  thin  annular  beam  in  a  long  drift  tube  is 

111  fe)  /  2/3  ,\3/2 


(35) 


WB 


2  In 


nl  (2/3  _ 


(36) 


a  result  that  prefectly  agrees  with  the  infinite-length  approximation  (see  Ref.  [3]).  The  finite- 
length  effects  are  shown  in  Fig.  1  which  compares  Eqs.  (32)  and  (35). 


V.3  Full  annular  beam 

The  final  limiting  case  is  a  full  annular  beam,  i.e.,  a  beam  completely  filling  the  drift  tube.  Under 
such  circumstances,  ri  =  r\  and  r0  =  Thus,  the  boundary  condition  coefficients  a\  and  b\  are 
zero.  Equation  (21)  becomes 

Jo  (fen)  No  {kr2)  -  Jo  (fcr2)  N0  {hr,)  =  0.  (37) 
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Figure  1:  Comparison  of  the  SCL  current  for  a  thin  annular  electron  beam  travelling  in  a  coaxial 
drift  tube  with  finite-length  effects  and  without.  The  outer  conductor  radius  is  0.01  cm. 


This  is  a  cross-product  of  circular  Bessel  functions  and  the  solution  k  can  be  approximated  (See 
Eqs.  [9.5.27]-[9.5.29]  of  Ref.  [4])  by 

7T  7~2  —  Tl 


k  S 

r2  —  r\  87rrir2 

Using  this  in  Eq.  (12),  the  SCL  current  upper  bound  for  a  full  annular  beam  is 


lUB 


2  2 
r2  —  rl 


7T 


r 2-n 


Lr2  -  ri  87rrir2J 


+ 


(7o2/3-l)3/2. 


(38) 


(39) 


VI  Concluding  Remarks 


A  theory  to  estimate  the  SCL  current,  including  finite-length  effects,  for  uniform,  relativistic 
electron  beams  propagating  in  circular  drfit  tubes  has  been  extended  to  include  annular  drift 
tubes  and  beams  of  arbitrary  radius  and  thickness.  For  limiting  beam  and  tube  geometries, 
closed-form  expressions  for  the  SCL  current  are  available.  For  other  cases,  a  numerical  solution 
must  be  used. 


8-11 


The  theory  assumes  a  strong  external  magnetic  field  to  prevent  radial  and  aximuthal  motion  of  the 
beam  particles.  Using  a  beam  of  uniform  density,  the  theory  is  not  self-consistent.  Furthermore, 
the  method  can  be  applied  only  to  those  cases  where  the  drift  tube  walls  are  grounded.  It  has 
recently  been  shown  that  biasing  the  inner  conductor  can  enhance  the  SCL  current,  see  Ref.  [3]. 
Methods  to  extend  the  present  theory  to  such  cases  is  underway. 
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